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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC) 
UNITS  OF  MEASUREMENT 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  (metric) 
units  as  follows: 


Multioh 

degrees 

feet 


0.01745329 

0.3048 


To  Obtain 
radians 
meters 
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SUMMARY 


This  report  describes  the  theory,  methodology,  and  verification  of  the  finite 
element  numerical  model  ADCIRC,  an  ADvanced  three-dimensional  CIRCulation  model 
developed  for  the  specific  purpose  of  generating  long  time  periods  of  hydrodynamic 
circulation  along  shelves,  coasts,  and  within  estuaries.  The  intent  of  the  model  is  to 
produce  long  numerical  simulations  (on  the  order  of  a  year)  for  very  large  computational 
domains  (for  example  the  entire  east  coast  of  the  US).  Therefore,  the  model  was 
designed  for  high  computational  efficiency  and  was  tested  extensively  for  both 
hydrodynamic  accuracy  and  numerical  stability.  The  results  of  these  tests  are  included 
in  this  report. 

The  ADCIRC  model  was  developed  by  the  Dredging  Research  Program  (a)  to 
provide  a  means  of  generating  a  database  of  harmonic  constituents  for  tidal  elevation 
and  current  at  discrete  locations  along  the  east,  west,  and  Gulf  of  Mexico  coasts,  and 
(b)  to  utilize  tropical  and  extratropical  global  boundary  conditions  to  compute  frequency 
indexed  storm  surge  hydrographs  along  the  US  coasts.  The  database  of  storm  and  tidal 
surface  elevation  and  current  data  is  being  developed  to  provide  site-specific 
hydrodynamic  boundary  conditions  for  use  in  analyzing  the  long-term  stability  of 
existing  or  proposed  dredge  material  disposal  sites. 

The  overall  intent  of  the  DRP  work  unit  is  to  provide  a  unified  and  systematic 
methodology  for  investigating  the  dispersive  or  nondispersive  characteristics  of  a  disposal 
site.  These  goals  can  be  realized  through  the  use  of  hydrodynamic,  sediment  transport, 
and  bathymetry  change  models.  The  ADCIRC  model  provides  the  tidal-  and  storm- 
related  hydrodynamic  forcings  necessary  for  site-specific  site  designation. 
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ADCIRC:  AN  ADVANCED  THREE-DIMENSIONAL  CIRCULATION  MODEL 


FOR  SHELVES.  COASTS.  AND  ESTUARIES 
THEORY  AND  METHODOLOGY 


PART  I;  INTRODUCTION 

1.  Interest  in  developing  a  more  accurate  technique  for  predicting  sea  surface 
elevation  and  circulation  in  coastal  areas  has  been  spurred  on  by  concerns  relating  to 
navigation,  shoreline  flooding,  pollutant  transport,  and  sediment  transport.  A  model 
for  computing  the  important  features  of  circulation  patterns  driven  by  tides,  wind, 
atmospheric  pressure  gradients,  and  ocean  currents  must  be  broad  in  scope  and  size. 

To  simplify  seaward  boundary  conditions,  yet  include  important  flow  details,  the  model 
must  encompass  large  domains  while  providing  a  high  degree  of  resolution  in  high- 
gradient  regions  as  well  as  in  nearshore  areas.  This  means  that  the  model  should 
allow  for  the  simultaneous  solution  of  flow  in  continental  shelf  regions,  coastal  areas, 
and  in  estuarine  systems.  The  model  should  solve  the  three-dimensional  conservation 
equations  [thereby  resolving  the  vertical  profile  of  horizontal  velocity]  instead  of  the 
widely  used  depth-integrated  conservation  equations.  This  is  necessary  since  it  is 
impossible  to  assume  a  relationship  between  bottom  stress  and  depth-averaged  velocity 
that  is  generally  valid  for  stratified  flows,  Ekman  layers,  and  wind-driven  circulation 
in  enclosed  or  semi-enclosed  basins  or  in  cases  where  wave  orbital  velocities  or 
suspended  sediment  concentration  gradients  are  significant  near  the  bottom. 

f  urthermore,  it  is  impossible  to  assume  values  for  momentum  dispersion  coefficients, 
which  are  inherent  in  depth-integrated  solutions,  that  are  generally  valid  in  complex 
flows. 

2.  The  requirements  of  very  large  domains,  a  high  degree  of  horizontal 
resolution  in  portions  of  the  domain,  and  the  resolution  of  i  ’.pidly  varying  vertical 
profiles  of  horizontal  velocity  place  strenuous  demands  on  even  the  largest 
supercomputers.  The  goal  in  the  development  of  ADCIRC  (ADvanced  three- 
dimensional  CIRCulation  model)  has  been  to  bring  together  algorithms  that  are  highly 
flexible,  accurate,  and  extremely  efficient.  These  issues  are  closely  interrelated  and 
have  been  emphasized  in  the  selection  of  discretization  techniques.  The  algorithms 
that  comprise  ADCIRC  allow  for  an  effective  minimization  in  the  required  number  of 
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degrees  of  freedom  for  a  desired  level  of  accuracy,  show  good  stability  characteristics, 
generate  no  spurious  artificial  modes,  have  minimal  inherent  artificial  numerical 
damping,  efficiently  separate  the  partial  differential  equations  into  small  systems  of 
algebraic  equations  with  time-independent  matrices,  and  are  capable  of  running  months 
to  years  of  simulation  while  providing  detailed  intra-tidal  computations. 

3.  The  framework  within  which  ADCIRC  has  been  developed  is  a  coupled 
external  mode  -  internal  mode  approach.  This  technique  has  proven  to  be  successful 
in  past  three-dimensional  models  and  can  significantly  reduce  the  cost  of  three- 
dimensional  hydrostatic  circulation  computations.  The  governing  equations  and  the 
basic  concept  behind  mode  splitting  are  discussed  in  detail  in  Part  II.  The  external 
mode  solution,  which  uses  the  well-known  depth-integrated  or  shallow-water  equations, 
is  discussed  in  Part  III.  Key  features  of  the  external  mode  solution  include  the  use  of 
a  generalized  wave-continuity  equation  (GWCE)  formulation  and  numerical 
discretizations  using  the  finite  element  (FE)  method  in  space  and  the  finite  difference 
(FD)  method  in  time.  Results  are  presented  using  the  external  mode  solution  as  a 
stand-alone,  two-dimensional  model  on  a  quarter  annular  test  case  and  the  North 
Sea/English  Channel  system.  Part  IV  focuses  on  the  internal  mode  solution'.  During 
the  development  of  ADCIRC,  a  novel  technique  was  discovered  that  replaces  velocity 
with  shear  stress  as  the  dependent  variable  in  the  internal  mode  equations.  The 
resulting  direct  stress  solution  [DSS]  allows  physically  realistic  boundary  layers  to  be 
included  explicitly  in  a  three-dimensional  model.  This  formulation  of  the  internal 
mode  equations  should  be  invaluable  for  modeling  coastal  and  shelf  circulation,  in 
which  the  bottom  and  surface  boundary  layers  comprise  a  significant  portion  of  the 
water  column,  and  for  modeling  processes  that  are  critically  dependent  on  boundary 
layer  physics  such  as  wave-current  interaction,  sediment  transport,  oil  spill  moveuient, 
ice  floe  movem  energy  dissipation,  physical-biological  couplings,  etc.  Thorough 
descriptions  of  the  DSS  formulation  and  testing  are  presented  in  Part  IV. 

4.  ADCIRC  is  being  developed  and  implemented  as  a  multi-level  hierarchy  of 
models.  A  2DDI  (two-dimensional,  depth-integrated)  option  solves  only  the  depth- 
integrated,  external  mode  equations  using  parametric  relationships  for  bottom  friction 
and  momentum  dispersion.  A  SDL  (three-dimensional,  local)  option  uses  horizontally 
decoupled  internal  mode  equations  to  solve  for  the  vertical  profile  of  horizontal 
velocity  and  to  evaluate  bottom  friction  and  momentum  dispersion  terms  for  the 
depth-integrated  external  mode  solution.  A  3DLB  (three-dimensional,  local, 
baroclinic)  option  includes  baroclinic  terms  as  a  diagnostic  feature.  Finally,  the  3D 
and  3DB  options  solve  the  complete  internal  mode  equations  for  nonstratified  and 


stratified  flows,  respectively.  At  present  ADCIRC-2DDI  is  fully  implemented  and 
operational,  ADCIRC-3DL  is  being  tested,  and  other  ADCIRC  versions  are  under 
development. 

5.  ADCIRC  achieves  a  high  level  of  simultaneous  regional/local  modeling, 
accuracy,  and  efficiency.  This  performance  is  a  consequence  of  the  extreme  grid 
flexibility,  the  optimized  governing  equation  formulations,  and  the  numerical  algorithms 
used  in  ADCIRC.  Together,  these  allow  ADCIRC  to  run  with  order  of  magnitude 
reductions  in  the  number  of  degrees  of  freedom  and  the  computational  costs  of  many 
presently  existing  circulation  models. 
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PART  II;  GOVERNING  EQUATIONS 


Three-dimensional  Equations  for  Nearly  Horizontal  Flow  in  Cartesian  Coordinates 


6.  A  survey  of  several  recent  review  volumes  (e.g.,  Heaps  1987;  Nihoul  and 
Jamart  1987;  ASCE  1988a, b;  Davies  1989)  indicates  that  the  turbulent  incompressible 
Reynolds  equations  simplified  using  the  Boussinesq  approximation  and  the  hydrostatic 
pressure  approximation  generally  form  the  basis  for  state-of-the-art  numerical  models 
of  coastal/shelf  circulation.  Although  these  equations  describe  fluid  motion  in  three 
dimensions,  because  of  the  simplification  of  the  vertical  momentum  equations,  they  are 
only  correct  for  nearly  horizontal  flow  (Koutitas  1987;  Abbot  1990).  Using  a 
right-handed  Cartesian  coordinate  system  these  equations  can  be  written  as 


du 


dvi 

W 


dv 


+  ^  +  It  “  ^ 


+  +  V. 


.dv 


-f  U^  +  V 


3y  ^ 
dv 


-dy 


,  dv 


(1) 

(2) 

(3) 


=  -  Pg 


(4) 


where 

f  =  2nsin^  =  Coriolis  parameter 
g  =  acceleration  of  gravity 
r  =  tide  generating  potential 
1/  =  molecular  viscosity 
p(x,y,z,t)  =  time-averaged  pressure 
/j(x,y,z,t)  =  density  of  water 
Po  =  reference  density  of  water 
t  =  time 


T  =  integration  time  scale  for  separating  turbulent  and  time-averaged  quantities 


1 

T 


1/  ^ 


1 

T 


i; 

r 


u'u'  dt  -  combined  viscous  and  turbulent  Reynolds  stress 
u'v'  dt  -  combined  viscous  and  turbulent  Reynolds  stress 
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T-zxC^.yi^.O  = 

r,y(x,y,z,t)  = 

7'yy(x,y,Z,t)  = 

'rzy(x.y»z.t)  = 


dw  1 
^  ■  T 


fT 

I 

0 

T 


u'w'  dt  -  combined  viscous  and  turbulent  Reynolds  stress 


1 

1  v'u'  dt  -  combined  viscous  and  turbulent  Reynolds  stress 

dv  1 

combined  viscous  and  turbulent  Reynolds  stress 

dw  1 

u  ~  T  1  “  combined  viscous  and  turbulent  Reynolds  stress 


(f>  =  degrees  latitude 

u(x,y,z,t),  v(x,y,z,t),  w(x,y,z,t)  =  time-averaged  velocities  in  the  x,  y  and  z  directions 

u'(x,y,z,t),  v'(x,y,z,t),  w'(x,y,z,t)  =  departures  of  the  instantaneous  turbulent 

velocities  from  the  time-averaged  velocities 

X,  y  =  horizontal  coordinate  directions 

z  =  vertical  coordinate  direction 

n  =  angular  speed  of  the  Earth  (7.29212xl0‘®  rad/s) 


7.  Using  the  vertical  momentum  equation,  pressure  can  be  eliminated  as  a 
dependent  variable  from  Equations  2  and  3,  to  give; 

du  ,  „5u  ,  „5u  .  ^ 


,  l/U  C/U  t 

^  +  v^  +  -  fv  - 


d 


-  r]  +  y^)  -  b.  +  m,  (6) 

^  +  ''1^  +  ^  I;  +  8^  ■  ^]  +  )  -  by  +  my  (6) 


where 


bx 


-  S_  ^ 


by  =  S-  I- 

^  Po  ^ 


(p-po)  dz  -  baroclinic  x  -  forcing 


(p-Po)  dz  -  baroclinic  y  -  forcing 


C(x,y>t)  =  free  surface  elevation  relative  to  the  geoid 


Mtx  =  T-, 
PoL 


^  ^ 

35c  by  j 


my  =  — 

^  Po 


-  horizontal  momentum  diffusion 

-  horizontal  momentum  diffusion 


Ps(x,y,t)  -  atmospheric  pressure  at  the  free  surface 


8.  The  solution  of  Equations  1,  5,  and  6  requires  the  following  boundary 
conditions: 
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a.  At  the  free  surface, 

Tzx  —  Tsx,  Tzy  —  Tsy  (J 

where  7’sx(x,y,t)  and  rsy(x,y,t)  are  wind  stresses  applied  at  the  water 
surface. 


b.  At  the  bottom, 

5h  .  „  5h  ,  5hl 

*  =  -  iK  +  “  af  +  ' 

rz%/Po  =  Tix/Po  =  kUb,  U,Ipo  =  'Tby/po  =  kVb 


(9) 

(10a) 


u  =  0,  V  =  0 


-h  -I-  Zo 


(10b) 


where  Tbx(x,y,t)  and  Tby(x,y,t)  are  bottom  stresses,  Ub(x,y,t)  and 
Vb(x,y,t)  are  near  bottom  velocities,  k  is  a  slip  coefficient  and  Zo  is 
the  effective  bottom  roughness  height  (e.g.,  Zq  =  kg/30  where  kg  is  the 
physical  bottom  roughness).  The  physic^y  correct  no-slip  condition, 
Equation  10b,  is  often  replaced  by  the  slip  condition.  Equation  10a,  to 
avoid  the  need  to  numerically  resolve  the  sharp  vertical  gradients  of  u 
and  V  that  exist  near  the  bottom.  A  quadratic  slip  conmtion  is 
obtained  by  setting 

k  =  Cd  (ug+  vg)*^^  (11) 

If  the  velocity  profile  is  logarithmic  between  the  elevation  where  Ub 
and  Vb  are  computed,  (-h+zb),  and  the  bottom,  (-h+Zo),  Cd  can  ^ 
defined  rigorously  as 

Cd  =  {i  ln[(zb-b)/(z,-h)l}-'  (12) 

where  k  is  the  von  Karman  constant.  Often  the  quadratic  slip 
condition  is  replaced  by  a  linear  slip  condition  by  setting  k  equal  to  a 
constant. 

c.  At  land  boundaries  normal  flux  is  specified.  Typically,  this  is  zero  for 
a  solid  boundary  or  nonzero  for  a  river  boundary. 

d.  At  open  boundaries  (either  along  the  ocean  or  at  rivers)  the  free 
surface  elevation,  C(x>yit)»  Is  specified,  a  radiation  boundary  condition 
is  used  to  allow  waves  to  enter  and  propagate  out  of  the  domain 
(Davies  and  Fumes  1980;  Rdd  1990),  or  the  discharge  is  specified. 


Three-dimensional  Equations  for  Nearly  Horizontal  Flow  in  cr  Coordinates 


9.  It  is  often  useful  to  transform  Equations  1,  5  and  6  into  a  bottom  and 
surface-following  "a"  coordinate  system.  By  this  means,  numerical  solutions  of  the 
transformed  equations  maintain  the  same  vertical  resolution  at  each  horizontal  grid 
point,  regardless  of  variations  in  depth  (Davies  1985;  Blumberg  and  Mellor  1987).  In 
a  general  <r-coordinate  system  (where  cr  =  a  at  the  free  surface  and  cr  =  b  at  the 
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bottom): 

=  X 

Va  =  y 

(T  =  a  + 
ta  =  t 
where 

H(x,y,t)  =  C  +  “  total  water  depth  to  the  free  surface 

h(x,y)  -  bathymetric  depth  relative  to  the  geoid 


(13a) 

(13b) 

(13c) 

(13d) 


(The  cr  subscript  is  used  to  denote  variables  iu  the  new  coordinate  system.) 


10.  Derivatives  are  converted  to  the  a-coordinate  system  using  the  chain  rule; 


d  _  d  ,  da  a  _  d  (a-b)  («<:  .  (<r-a)  8H  1  d 

m  -  -d^w-  m-  ^ 


(14a) 


d  _  d  ,  da  d  (a— b)  ,  (<y-a)  5H  1  3 

^  W  H  (a-b)  ^ 


(14b) 


d  (a-b)  d 


(14c) 


d  _  d  ,  da  d  _  d  (a— b)  f^C  ■  (^^-a)  5H  1  5 

M  =  m  W  =  'dt~  Tn  ^ 


(14d) 


11.  The  velocity  component  aligned  in  the  a  direction  is  defined  as 


do  ^  (a^) 
""  -  3t.  TT 


(a-b)  3r, 


I  ■  -B.*  IS)  g.) 


12.  The  baroclinic  forcings  bx<,,  by„,  and  the  horizontal  momentum  diffusion 
mxg,  my^  in  the  a  coordinate  system  become: 


bxa  -  2^^  %+  (p-Po)  dor]  +  (cr-a)  (16a) 

L  a 

ha  ==  (p-Po)  +  (o^-a)  ^J^P-Po)  (16b) 

I  <F 
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13.  Substituting  Equations  13  -  17  into  Equations  1,  and  4-6,  and 
rearranging  terms  gives  the  three-dimensional  governing  equations  in  the  a-coordinate 
system.  Dropping  the  a  subscripts  for  notational  convenience,  the  transformed 
equations  are 


5uH  ,  dvE  , 

(18) 

dn  .  du  , 

1 

II 

+  ^  -  b.  + 

(19) 

dv 

+ 

(20) 

(21) 

14.  The  a  equations  use  the  same  boundary  conditions  as  the  original 
equations  with  the  exception  that  w^=  0  at  the  free  surface  and  at  the  bottom. 


Vertically  Integrated.  Two-dimensional  Equations  for  Nearly  Horizontal  Flow 


15.  The  three-dimensional  equations  can  be  integrated  over  the  vertical  to 
yield  a  set  of  two-dimensional  equations  for  free  surface  displacement  and 
depth-averaged  velocity.  In  conservative  form  these  equations  are; 


aJH 


+ 


5VH 

W 


(22) 


auH  ,  auuH  .  5UVH 

w~  +  ~mr  ~w~ 


-  fVH  =  -  H 


+  g(C  -  a,) 
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(23) 


+  M.  +  D.  +  B.  +  ISi  -  32 


Po  Pq 


m 


5y 


oy[po 

+  My  +  Dy  +  By  +  ^  ^ 

^  ^  ^  Po  Po 


(24) 


where 


a  =  eflective  Earth  elasticity  factor  (a  =  0.69) 

Bx  s  -  I  bx  dz  -  depth^ntegrated  baroclinic  forcing 


-h 


B 


by  dz  -  depth-integrated  baroclinic  forcing 
-h 

~  ~  ~  momentum  dispersion 


5Duv  dPw 


-  momentum  dispersion 


Duu  =  [  uu  dz,  Duv  =  [  vu  dz,  Dw  =  [  vv 
-h  -h  -h 

r;(x,y,t)  -  Newtonian  equilibrium  tide  potential 


dz 


Mx  = 


^  f 


-h  “h 


dz  -  depth-integrated,  horizontal  momentum  diffusion 


n  Q 

My  =  ^  dz  +  ^  dz  -  depth-integrated,  horizontal  momentum  diffusion 

-h  ^ -h 

1 

U(x,y,t)  =  w  u  dz  -  depth-averaged  horizontal  velocity 
-h 

1 

V(x,y,t)  =  ^  V  dz  -  depth-averaged  horizontal  velocity 


u(x,y,z,t)  =  u  -  U  -  departure  of  horizontal  velocity  from  depth-averaged  velocity 
v(x,y,z,t)  =  V  -  V  -  departure  of  horizontal  velocity  from  depth-averaged  velocity 

16.  In  non-conservative  form,  the  vertically  integrated  momentum  conservation 
equations  are: 

+  u|2  +  vH  -  fV  =  -  +  g«  -  a,)] 

+  h[««  +  +  b.  +  ^  -  a.]  (25) 
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(26) 


It  +  +  n’  =  - +  g(C  -  0^)] 

+  ^[m,  +  d,  +  b,  +  -  ^] 

17.  The  derivation  of  the  Newtonian  equilibrium  tide  potential,  t),  is  presented 
by  Reid  (1990).  A  practical  expression  for  rj  given  by  Reid  is 

1]{X,(p,t)  =  J  Cjn  fjn(to)  hj(^)  COS[2jr(t-to)/Tjn  +  jA  +  ^n(to)]  (27) 

n.j 

where 

Cjn  =  constant  characterizing  the  amplitude  of  constituent  n  of  species  j  (Table  1) 
fjn(t)  =  time-dependent  nodal  factor 

j  =  0,  1,  2  -  tidal  species  (j=0  declinational,  j=l  diurnal,  j=2  semidiurnal) 

Lo  =  3  sin2(^)  -  1 
Li  =  sin(2^) 

L2  =  cos  2(0) 

A,  0  =  degrees  of  longitude  and  latitude,  respectively 
to  =  reference  time 

Tjn  =  constaiit  characterizing  the  period  of  constituent  n  of  species  j  (Table  1) 

^n(t)  -  time-dependent  astronomical  argument 

Values  for  fjn  and  can  be  computed  from  tables  (e.g.,  Schureman  1941)  or  using 
available  harmonic  analysis  packages  (e.g.,  Foreman  1977). 

18.  The  gradient  of  ar)  results  in  the  effective  tide-producing  force.  The 
factor  a  accounts  for  the  reduction  in  the  field  of  gravity  due  to  the  existence  of 
small  tidal  deformations  of  the  Earth’s  surface  called  Earth  tides.  The  value 

a  =  0.69  is  the  ratio  of  the  theoretical  period  of  the  Earth’s  wobble  derived  by  Euler 
(assuming  the  Earth  to  be  a  perfectly  rigid  sphere)  to  the  observed  period  of  the 
Earth’s  wobble  (Reid  1990).  (Therefore  a  is  a  global  measure  of  the  rigidity  of  the 
Earth.  For  reference,  a  =  1  would  correspond  to  a  perfectly  rigid  sphere.)  a  =  0.69 
has  been  used  for  modeling  global  ocean  tides  by  investigators  including  Schwiderski 
(1980)  and  Hendershott  (1981). 

19.  Due  to  their  computational  efficiency,  models  based  on  the  vertically 
integrated  equations  have  been  widely  used  for  modeling  coastal,  shelf,  and  even  open 
ocean  circulation  (e.g.,  Leendertse  1967;  Wang  and  Connor  1975;  Spaulding  1984; 

Smith  and  Cheng  1987;  Werner  and  Lynch  1987;  Walters  1987;  Vincent  and  Le 
Provost  1988;  Westerink,  Stolzenbach,  and  Connor  1989;  Signell  1989).  All  of  the 
physics  contained  in  the  original  three-dimensional  governing  equations  are  embedded 
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Table  1 

Constants  for  the  Prindpal  Tidal  Constituents  Cfrom  Reid  19901 


C 

T 

Species 

Constituent 

m 

solar  days  or  hrs* ** 

0 

Mf 

fortnightly  lunar 

0.041742 

13.660791d 

M„ 

monthly  lunar 

0.022026 

27.554553d 

Ssa 

semiannual  solar 

0.019446 

182.6211d 

Sa 

annual  solar 

♦* 

365.2597d 

1 

Ki 

luni-solar 

0.141565 

23.9344696h 

0, 

principal  lunar 

0.100514 

25.8193417h 

Pi 

principal  solar 

0.046843 

24.0658902h 

Qi 

elliptical  lunar 

0.019256 

26.8683566h 

2 

Mj 

principal  lunar 

0.242334 

12.4206012h 

S2 

principal  solar 

0.112841 

12.0000000h 

N2 

elliptical  lunar 

0.046398 

12.6583482h 

K2 

luni-solar 

0.030704 

11.9672348h 

*One  lunar  day  =  1.035050  solar  days  or  24.8412  solar  hours 

**The  annual  solar  tide  is  heavily  dependent  on  seasonal  heating  and  cooling  of  the 
ocean,  as  well  as  radiation  pressure. 
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in  the  vertically  integrated  equations  if  the  bottom  stress  and  the  momentum 
dispersion  terms  are  specified  correctly.  Although  more  sophisticated  approaches  have 
been  developed  for  specialized  conditions  (Lynch  and  Officer  1985;  Nihoul  and  Djenidi 
1987;  Tee  198’5’;  Poon  1988;  Jenter  and  Madsen  1989),  bottom  stress  is  usually 
parameterized  as  a  collinear  function  of  the  depth-averaged  velocity,  and  momentum 
dispersion  is  either  neglected  or  represented  as  a  "diffusion-like''  function  of  the 
depth-averaged  velocity  (Bedford  1984). 

20.  Parameterized  bottom  stress  relationships  are  typically  quadratic  in  the 


depth-averaged  velocity  and  of  the  form 

^  =  Cf  (U2  +  U 

Po 

(28a) 

^  =  Cf  (U2  +  V2)i/2  Y 

Po  '  ' 

(28b) 

where  Cf  is  computed  using  one  of  the  following  relationships: 

II 

O 

(29a) 

C’ 

(29b) 

C,  = 

(29c) 

In  Equation  29,  is  the  Darcy-Weisbach  friction  factor,  C  is  the  Chezy  friction 

coefficient,  and  n  is  the  Manning  friction  factor. 

21.  The  depth-integrated  lateral  momentum  diffusion  terms  are  typically 
lumped  together  with  the  momentum  dispersion  terms  into  a  standard  isotropic  and 
homogeneous  eddy  diffusion/dispersion  model  (Blumberg  and  Mellor  1987) 


Mx  +  Dx  = 
My  +  Dy  = 


d^m 


+ 


dy  2  dx.  d  y 


pMD 

hh  j 


ra^vH 

dx  2 


+  2 


d^VH 
dy  2 


a^UHl 


(30a) 

(30b) 


where  is  a  horizontal  eddy  diffusion/dispersion  coefficient.  Equation  30  is  based 

directly  on  a  molecular  diffusion  analogy  as  applied  to  depth-integrated  flow.  Kolar 
and  Gray  (1990)  use  a  slightly  simpler  model  that  approximates  Equation  30  as: 


Mx  +  Dx 


r^^uH 


dy  2 


(31a) 
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My  +  Dy  = 


t^MD 

i^h2 


r^^vH 

~ur^ 


a^VHl 


(31b) 


where  e“  j 


is  an  eddy  diffusion/dispersion  coefficient  that  will  generally  not  be  equal 


to  E 


MD 

hr 


22.  For  flows  with  horizontal  length  scales  that  are  large  compared  to  the 
depth,  Mx  and  My  are  negligible  in  the  momentum  balance  in  Equations  25  and  26 
(Blumberg  and  Mellor  1987).  D*  and  Dy  are  similarly  small  when  the  velocity  profile 
is  nearly  uniform  over  the  vertical.  In  such  flows  Eh^  or  Ehj  are  either  set  to  zero 
or  kept  at  a  relatively  small  value  to  provide  stability  to  the  numerical  scheme.  (The 
latter  must  be  done  with  considerable  caution  to  ensure  that  the  contributions  of  these 
terms  in  the  momentum  equations  remain  small.  Otherwise,  the  model  solutions  will 
be  artificially  altered.)  Conversely,  when  the  velocity  profile  varies  strongly  over  the 
vertical,  Dx  and  Dy  may  have  a  significant  contribution  to  the  momentum  balance. 

23.  For  tidad  flows  in  relatively  shallow,  unstratified  waters,  depth-4ntegrated 
computations  that  make  use  of  the  parameterizations  given  in  Equations  28-31 
appear  to  work  reasonably  well  (although  detailed  studies  of  tidal  constituent  dynamics 
indicate  that  all  of  the  flow  physics  are  not  captured  in  two-dimensional  simulations 
due  to  the  form  of  the  bottom  friction  term  (Westerink,  Stolzenbach,  and 

Connor  1989)).  However,  in  wind-driven  flows,  stratified  flows,  Ekman  layers,  or 
when  wave  orbital  velocities  or  suspended  sediment  gradients  are  significant  near  the 
bottom,  the  simple  parameterizations  for  bottom  friction  and  momentum  dispersion 
given  »above  become  entirely  inadequate.  Also,  since  the  depth-averaged  velocity  may 
be  very  different  from  the  actual  velocity  at  a  specific  elevation  in  th?  water  column 
(particularly  if  flow  reversal  occurs  over  the  depth),  the  use  of  the  depth-averaged 
velocity  in  a  transport  model  (e.g.,  for  sediment  transport)  may  cause  considerable 
error  in  predicted  transport  patterns.  Therefore,  for  many  applications  of  practical 
interest,  a  model  based  solely  on  the  vertically  integrated  governing  equations  is  not 
adequate. 


Mode  Splitting 

24.  Unfortunately,  numerical  solutions  of  the  three-dimensional  governing 
equations  require  substantially  increased  computer  time  and  storage  in  comparison  to 
solutions  of  the  vertically  integrated  equations.  To  help  minimize  this  cost,  most 
three-dimensional  models  use  some  type  of  mode-splitting  scheme.  Mode  splitting  is 
accomplished  by  solving  the  two-dimensional,  vertically  integrated,  "external  mode" 
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equations  for  the  free  surface  displacement  (and  sometimes  the  depth-averaged 
velocity).  The  external  mode  solutions  are  then  used  to  force  "internal  mode" 
equations  that  account  for  the  vertical  propagation  of  momentum.  The  internal  mode 
equations  are  solved  for  the  vertical  profile  of  velocity  and  the  results  used  to 
compute  Tbx,  Tby,  Dx,  and  Dy  for  subsequent  external  mode  calculations.  Internal 
mode  equations  have  been  generated  by  integrating  the  three-dimensional  equations 
over  discrete  layers  in  the  vertical  and  then  subtracting  the  equations  for  adjacent 
layers  (Simmons  1974;  Sheng  and  Lick  1980),  by  subtracting  the  external  mode 
equations  from  the  three-dimensional  equations  (Wang  1982;  Sheng  1983;  Davies  1985), 
by  differentiating  the  three-dimensional  equations  in  the  vertical  direction  (Tee  1979), 
or  by  using  the  three-dimensional  equations  themselves  (Blum berg  and  Mellor  1987; 
Lynch  and  Werner  1991).  (The  internal  mode  equations  and  their  solution  are 
discussed  in  detail  in  Part  IV  of  this  report.)  Mode  splitting  allows  the  free  surface 
elevation  to  be  evaluated  with  the  computational  efficiency  of  a  vertically  integrated 
model.  This  can  be  quite  important  since  the  allowable  time  step  for  this 
computation  is  often  severely  constrained  by  accuracy  requirements  or  a  Courant 
stability  criterion.  Since  the  internal  mode  calculations  are  free  from  surface  gravity 
waves,  the  vertical  profile  of  velocity  can  often  be  computed  using  a  significantly 
larger  time  step  than  the  free  surface  elevation. 

25.  In  effect,  mode  splitting  replaces  the  parameterizations  of  bottom  stress 
and  momentum  dispersion  used  in  a  purely  two-dimensional  model  with  values 
computed  from  the  vertical  profiles  of  velocity  generated  by  the  internal  mode 
equations.  Therefore,  the  vertically  integrated,  external  mode  computations  do  not 
require  parameterizations  of  either  bottom  stress  or  momentum  dispersion  in  terms  of 
the  depth-averaged  velocity.  The  only  parameterizations  maintained  in  the  external 
mode  equations  are  for  the  horizontal  momentum  diffusion  terms.  These  terms  are 
usually  insignificant  in  the  momentum  balance,  although  for  small-scale  computations 
horizontal  momentum  diffusion  can  be  a  physically  important  process.  Most  often  the 
horizontal  momentum  diffusion  terms  are  retained  only  to  provide  numerical  stability 
and  are  parameterized  with  expressions  identical  to  Equations  30  and  31,  i.e.. 
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or  alternatively. 
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(33b) 


My  =  E^2 

where  E^j 


ra^vH  .  a^VHl 

and  E^j  are  eddy  coefficients  for  horizontal  momentum  diffusion. 


Vertical  Turbulent  Closure 

26.  The  internal  mode  equations  require  the  parameterization  of  the  vertical 
turbulent  momentum  transport  terms,  and  ray,  (also  called  the  vertical  shear 
stresses).  These  terms  can  dominate  the  momentum  balance  in  portions  of  the 
domain  and  it  is  therefore  critical  to  find  an  adequate  closure  scheme.  Turbulent 
closure  has  been  and  continues  to  be  the  subject  of  considerable  research.  Recent 
summaries  of  this  work  include  Mellor  and  Yamada  (1982);  Rodi  (1984,  1987); 

Ferziger  (1987);  Johns  and  Oguz  (1987);  and  ASCE  ( 1988a, b).  The  most  general 
approach  is  to  solve  transport  equations  for  the  turbulent  velocity  correlations  that 
make  up  the  turbulent  stresses  (stress/flux  models).  However,  this  adds  considerably 
to  the  computational  burden  of  a  three-dimensional  model.  Models  based  on  this 
technique  have  had  little  testing  and  virtually  no  application  to  geophysical  flows 
(ASCE  1988b).  Also,  it  appears  that  these  models  offer  no  decisive  advantage  in 
shear  flows  (Launder  1984).  Alternatively,  the  vertical  shear  stresses  can  be 
parameterized  in  terms  of  the  mean  velocity  field  using  eddy  viscosity  relationships  of 


the  form 

Tzx  _  p  ^ 

(34a) 

=  Ey 

po  ^  m 

(34b) 

On  dimensional  grounds  the  vertical  eddy  viscosity  Ey  should  be  proportional  to  a 
velocity  scale  v  multiplied  by  a  length  scale  I,  both  of  which  are  characteristic  of  the 
turbulent  motion.  Particularly  simple  expressions  such  as  the  Prandtl  mixing  length 
model  can  be  found  for  v  and  I  for  boundary-layer  type  flows  (Rodi  1987).  In  more 
complex  flows,  v  has  been  related  to  the  square  root  of  the  total  turbulent  kinetic 
energy,  k.  The  terms  k  and  /  (or  some  combination  of  k  and  I  such  as  can 

be  solved  for  using  quasi-empirical  transport  equations  or  specified  using  empirical 
algebraic  expressions.  The  primary  limitations  to  the  eddy  viscosity  approach  are  its 
inability  to  simulate  counter  gradient  transport  or  to  a^'count  for  nonisotropic 
turbulence.  A  third  choice  for  expressing  the  turbulent  stresses  lies  between  the 
stress/flux  models  and  the  eddy  viscosity  models  in  complexity  and  potential  for 
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representing  co>^.Dlex  flows.  In  this  approach  algebraic  expressions  (approximations  to 
the  transport  equations  used  in  stress/flux  models)  relate  the  vertical  stresses  to  k  and 
I  (or  e)  without  the  use  of  an  eddy  viscosity  hypothesis. 

27.  Eddy  viscosity  models  are  by  far  the  most  widely  used  method  for 
representing  vertical  momentum  transport  in  coastal  flows.  These  models  can  be 
expected  to  work  reasonably  well  in  such  applications,  since  the  water  column  is 
typically  dominated  by  the  bottom  and  surface  boundary  layers. 
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PART  III:  EXTERNAL  MODE  SOLUTION 


Selection  Considerations  for  the  External  Mode  Solntion 

28.  A  basic  objective  in  the  development  of  ADCIRC  is  to  provide  the  abibty 
to  perform  computations  on  very  large  domains.  This  requires  selecting  algorithms 
that  satisfy  interrelated  requirements  of  a  high  level  of  grid  flexibility,  accuracy,  and 
efficiency.  To  ensure  a  high  degree  of  solution  accuracy,  the  discretization  scheme 
must  have  numerical  amplitude  and  phase  propagation  characteristics  that  are  nearly 
identical  to  the  analytical  characteristics  even  for  relatively  poorly  resolved 
wavelengths  (e.g.,  good  correspondence  down  to  at  least  A/ Ax  =  20,  where  A  is  the 
wavelength  and  Ax  the  grid  spacing).  Furthermore,  solution  accuracy  requires  that  all 
wavelengths  with  significant  energy,  (e.g.,  as  generated  in  regions  of  rapidly  varying 
flow,  geometry,  and/or  topography),  be  well-resolved.  A  high  degree  of  solution 
efficiency  requires  that  the  algorithm  minimizes  both  the  number  of  degrees  of  freedom 
and  the  operations  required  per  degree  of  fireedom  per  time  step.  Minimization  of  the 
number  of  degrees  of  freedom  is  constrained  by  the  need  to  provide  resolution  on  a 
localized  basis  and  is  highly  dependent  on  the  accuracy  and  the  grid  flexibibty  of  the 
numerical  scheme. 

29.  Because  grid  flexibility  is  pivotal  to  solution  accuracy  and  efficiency, 
various  strategies  have  been  devised  to  allow  variations  in  grid  size  over  a  model 
domain,  A  nested  grid  approach  offers  one  solution.  However,  unless  the  grids  are 
coupled,  this  approach  cannot  properly  account  for  flow  interactions  between  the 
various  grids.  Stretched  FD  grids  offer  the  possibility  of  providing  local  refinement 
within  a  single  grid.  However,  cell  aspect  ratio  requirements  limit  the  degree  of  grid 
size  variability.  Furthermore,  since  cell  size  in  the  x  direction  is  fixed  for  all  y 
locations  for  a  given  x  and  vice  versa,  portions  of  the  domain  are  often  over-refined. 
Boundary-fitted  FD  schemes  that  utilize  conformal  mapping  techniques  allow  the  land 
boundaries  to  be  well-represented  in  addition  to  offering  local  refinement  possibilities. 
However,  these  techniques  suffer  from  the  same  shortcomings  as  stretched  FD 
approximations  and  often  significant  difficulties  are  encountered  in  finding  a  suitable 
transformation  function  for  complex  geographic  regions.  The  FE  algorithms  based  on 
triangular  elements  are  highly  flexible  and  can  provide  local  refinement  in  a  systematic 
and  optimal  fashion.  In  fact,  circulation  computations  for  tides  and  storm  surge  in 
the  Gulf  of  Mexico  (Westerink  et  al.,  in  press)  have  been  achieved  with  cell  area 
ratios  greater  than  1  to  15,000. 
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30.  Algorithm  accuracy  per  degree  of  freedom  is  another  critical  issue  in  the 
selection  of  an  external  mode  solution  algorithm.  FD  schemes  were  successful  fairly 
early  in  their  development  owing  to  the  use  of  the  staggered  or  C  grid  approach 
(Hansen  1956;  Leendertse  1967).  Early  FE  schemes  were  plagued  with  severe  spurious 
modes  that  required  the  heavy-handed  addition  of  non-physical  dissipation  and 
resulted  in  very  poor  accuracy  characteristics  (Gray  1982).  It  was  not  until  the 
introduction  of  the  wave-continuity  equation  (WCE)  formulation  that  robust  and 
highly  accurate  FE  schemes  emerged  (Lynch  and  Gray  1979).  The  WCE  formulation 
is  based  on  the  rearrangement  of  the  continuum  equations  prior  to  any  spatial 
discretization.  Extensive  numerical  testing  has  demonstrated  that  FE-based  WCE 
solutions  produce  very  accurate  results  (Lynch  and  Gray  1979;  Lynch  1981;  Walters 
and  Carey  1983;  Walters  1983  and  1984).  It  has  also  been  shown  that  the 
fundamental  success  of  the  WCE  FE  scheme  lies  in  its  ability  to  propagate  2 Ax 
waves  (Platzman  1981;  Foreman  1983).  (This  is  also  the  reason  why  the  C  grid  FD 
solutions  are  successful.) 

31.  Finally,  overall  algorithm  efficiency  is  essential  in  the  selection  of  an 
external  mode  solution.  In  general,  implicit  methods  are  more  useful  in  long  wave 
computations  than  explicit  methods,  particularly  when  small  cells  or  elements  are  used. 
However,  the  use  of  implicit  methods  typically  results  in  time-dependent  matrices  that 
must  be  reassembled  and  re-solved  at  every  time  step.  This  increases  the 
computational  burden  significantly.  The  FD  methods  overcome  this  problem  by 
implementing  an  alternating  direction  implicit  (ADI)  type  approach  that  reduces  a 
two-dimensional  problem  to  a  sequence  of  one-dimensional  problems,  resulting  in 
significant  computational  savings  for  large  problems.  It  is  not  possible  to  apply  the 
ADI  approach  to  FE-based  methods.  However,  a  WCE  FE-based  solution  has  been 
formulated  that  decouples  the  solutions  for  elevation  and  velocity  and  allows  the  use 
of  time-independent  matrices  for  the  elevation  solution  and  diagonal  matrices  for  the 
velocity  solution.  These  features  have  produced  a  highly  efficient  WCE  FE  solution 
called  the  generalized  wave-continuity  equation  (GWCE)  formulation  (Kinmark  1985). 

32.  Careful  consideration  of  the  requirements  for  grid  flexibility  and  a  high 
level  of  accuracy  and  efficiency  led  to  the  selection  of  the  FE-based  GWCE 
formulation  for  the  external  mode  solution  in  ADCIRC.  Extensive  analysis,  testing, 
and  field  applications  of  the  GWCE  during  the  past  decade  have  demonstrated  the 
unparalleled  capabilities  and  robustness  of  the  scheme. 
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Development  of  the  Generalized  Wave-Continuitv  Equation 


33.  The  GWCE  formulation  is  a  specifically  designed  WCE  formulation  that 
yields  a  discrete  system  of  equations  with  time-nndependent  matrices.  Time- 
independent  system  matrices  are  critical  in  minimizing  the  computational  cost  for 
finite-element-based  solutions  due  to  the  expense  of  both  the  matrix  assembly  and 
decomposition  steps.  The  GWCE  is  based  on  the  primitive  depth-integrated 
continuity  equation,  Equation  22,  and  the  primitive  depth-integrated  conservation  of 
momentum  equations  in  conservative  form.  Equations  23  and  24.  The  primitive 
continuity  equation  is  differentiated  vrith  respect  to  time  to  yield: 


^  ^t^  wdj  -  ^ 


The  primitive  momentum  equations  are  differentiated  with  respect  to  x  and  y, 
respectively,  and  rearranged  as: 
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Equations  36  and  37  are  then  substituted  into  Equation  35: 
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Finally,  the  primitive  continuity  equation  is  multiplied  by  a  constant,  Tq,  and  added 
to  Equation  38: 
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-  H  ^  1^  +  g«  -  <«,))  +  M,  +  D,  +  By  +  +  ToVH}  =  0  (39) 

34.  The  advective  terms  in  Equation  39  are  in  conservative  form.  Our 

experience  indicates  that  if  these  terms  are  put  into  non-conservative  form,  improved 
numerical  stability  is  obtained  when  advection  is  dominant  in  the  global  or  local  force 

balance.  The  advective  terms  in  the  GWCE  are  reformulated  by  expanding  the 

derivatives  and  substituting  in  the  primitive  continuity  equation,  Equation  22. 
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35.  The  lateral  closure  model  in  ADCIRC  is  the  simplified  eddy  viscosity 
model  of  Kolar  and  Gray  (1990),  Equation  33.  Substituting  this  into  Equation  40 
gives: 
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where  Ej^j  is  the  generalized  lateral  diffusion/dispersion  coefficient.  For  the  2DDI 

option,  E}j2  represents  the  combined  effects  of  both  lateral  diffusion  and  dispersion. 

Therefore  Ejjj  =  and  D*  and  Dy  are  both  set  to  zero.  For  the  three-dimensional 

ADCIRC  options,  Ej^j  represents  only  lateral  diffusion.  In  these  cases,  Ejjj  =  E^j, 

and  Dx  and  Dy  are  explicitly  computed  from  the  internal  mode  solution.  It  is 
assumed  that  E^^j  is  constant  in  time  and  space  and  that  it  has  a  value  of  zero  on 

the  boundaries  of  the  domain. 

36.  The  lateral  diffusive/dispersive  terms  in  Equation  41  can  be  conveniently 
rearranged  to  decrease  the  functional  continuity  requirements  for  the  symmetrical  weak 
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weighted  residual  formulation  from  back  to  C®  as  is  the  case  for  the  GWCE 
formulation  without  any  lateral  closure  model  (Kolar  and  Gray  1990).  Rearranging 
the  spatial  derivatives  of  the  lateral  diffusive/dispersive  terms  in  Equation  41  gives: 

0  It  ■  H  ~  I7  +  ®  H 

+  D,  +  B.  +  -h  roUH}  +  I  {V  1^  -  UH  1^  -  VH  1^  -  fUH 
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+  Ehj  (^3^  +  ^)]  +  Ehj  +  ^)]  =  0  (42) 

The  primitive  continuity  equation,  Equation  22,  can  be  used  to  substitute  for  the 
divergence  of  flux  in  the  lateral  diffusion/dispersion  terms  in  Equation  42  to  give: 
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37.  Equation  43  can  be  solved  in  conjunction  with  the  primitive  conservation 
of  momentum  equations  in  either  conservative  or  non-conservative  form.  ADCIRC 
uses  the  non-conservative  momentum  equations,  Equations  25  and  26.  Incorporating 
the  same  simplified  eddy  viscosity  model  into  the  non-conservative  momentum 
equations  gives: 
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38.  The  bottom  stress  in  Equations  43  -  45  is  expressed  using  a  drag  tensor 
similar  to  that  proposed  by  Jenter  and  Madsen  (1989): 


and  7  is  the  angle  measured  counter  clockwise  &om  the  depth-averaged  velocity  vector 
to  the  bottom  stress  vector. 

39.  Defining 

P=  f  +  r*sin(7)  (47a) 

ri  =  r*cos(7)  (47b) 

aiiH  substituting  Equations  46  and  47  into  Equations  43  —  45  gives  the  GWCE  and 
momentum  equations  in  final  form: 
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40.  In  the  2DDI  option,  the  bottom  stress  and  depth-averaged  velocity  are 
assumed  to  be  co-4inear  (7  =  0).  Cf  is  specified  directly  as  an  input  parameter  or 
computed  using  one  of  the  relationships  given  in  Equation  29.  In  the  three- 


dimensional  ADCIRC  options,  7  and  Cf  are  computed  using  rb*  aJ^d  rby  from  the 
internal  mode  solution.  As  noted  above,  7  is  the  angle  measured  counterclockwise 
from  the  depth-averaged  velocity  vector  to  the  bottom  stress  vector.  Cf  is  determined 


as: 


Po(U‘  +  V’) 


(51) 


It  is  easily  shown  that  Equations  46,  47,  and  51  introduce  the  bottom  stresses 
computed  in  the  internal  mode  solution  directly  into  the  external  mode  equations. 


Development  of  Weighted  Residual  Statement 


41.  To  develop  a  Galerkin  weighted  residual  statement  for  the  GWCE,  Equation 
48  is  weighted  by  the  interpolating  basis  function,  and  spatially  integrated  over 
the  interior  domain,  D,  giving: 

^  -N  (52) 

where 

<o,6>n  =  ^/  0  6  dn 
n  =  the  global  domain 

N  =  number  of  nodes  in  the  spatial  discretization 

A.  =  U  If  -  UH  |2  -  VH  ^  +  f'VH  -  H  (^  +  g«  -  a,)] 

-  Ek,  1^  +  D,  +  B,  +  +  (r<,-Ti)UH  (53a) 

A,  =  V  If  -  UH  1^  -  VH  f  -  f'UH  -  H  (Ei  +  g(C  -  a,)] 

-  Ehl  +  ^  +  (ro-rt)VH  (53b) 

Applying  Gauss’s  theorem  to  the  integrals  in  Equation  52  that  contain  spatial 
derivatives  gives: 

<|t^>  ~  -  <Ay,  |^>n 

=  -  j,  [Ax«nx  +  Mnyl^i  "  N  (54) 
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where  F  is  the  boundary  of  the  domain  fi.  The  direction  cosines  are  defined  as: 


Onx  =  cos{^x)  (55a) 

Ofay  =  cos(^y)  (55b) 


where  and  &y  are  the  (spatially  varying)  angles  measured  to  the  outward  normal  at 
any  point  along  the  boundary  from  the  positive  x  and  y  axes,  respectively. 

42.  Using  the  conservative  form  of  the  momentum  equations,  Equations  23  and 
24,  recasting  the  advective  terms  in  Equation  53  into  conservative  form,  and  using  the 
simplified  lateral  diffusion  model,  Equation  33,  Ax  and  Ay  can  be  written  as: 


A,  = 


5UH 


-  ^h2  + 


Ay  -  - l^ha  + 


d^VH. 

'dp-^  ~ 


+  ^oUH 
+  r,VE 


(56a) 

(56b) 


Substituting  Equation  56  into  the  line  integral  in  Equation  54  and  assuming  that  Ejjj 
is  zero  on  the  boundary.  Equation  54  becomes: 

It  “  lx  ~  ly^ 

“  ®  S  8(C  -  <w?)l  -  Et,  ^  +  D,  +  Bx  +  ^>n 

-  <V  ^  -  UH  -  VH  f  -  f'UH  -  H  I;  (fe  +  6«  -  a,)]  -  E„  ^ 

+  Dy  +  By  +  +  (r<,-ri)VH.  ^>a  =  -  I  ||{UHa„x  +  VHa„,) 

+  To(UHOfnx  +  VHttny)]^!  dF  i  =  1,  ...N  (57) 

43.  The  terms  that  involve  partial  derivatives  of  the  barometric  pressure, 
surface  elevation,  and  Newtonian  equilibrium  tidal  potential  can  be  written  as: 

“  If  +  8(<  -  Of)]  =  Sk  ^  +  I  ^  +  6H  (^  -  OT,)  (58a) 

H  I  -  “^)1  =  Si  I  +  §  I  +  SH  I  (^  -  «5)  (58b) 

44.  The  normal  flux  across  the  boundary  is  defined  as: 

Qn  =  UHOnx  "b  VHOfny  (59) 

45.  The  line  integral  in  Equation  57  is  non-zero  only  on  flux-specified 
boundaries,  F^.  Using  the  specified  normal  flux  Qn*  for  Qn,  and  substituting 

Equations  58  and  59  into  Equation  57  gives  the  final  symmetrical  weak  weighted 


31 


residual  statement  for  the  GWCE: 


+  <To^,<t>i>n  +  <g^^>  ^>n  +  ^>n  +  ^>n 

+  Eha<^,  ^>n  =  <u|^.  ^>n  +  <v|^.  ^>n  +  <Wx,  ^>n 
+  <Wy,  ^>n  -  ^  (#"*  +  ^oQn*)  dr  i  =  1,  ...N  I 


where 


w,  ^  -  UH  H  -  VH  1^  +  f' VH  -  I  g!  -  gH  (-^  -  a,)  +  D,  +  B. 


+  ^  +  (ro-ri)UH 


(6U) 


W,  .  -  UH  W  -  VH  1^  -  f'UH  -  §  gH  -  an)  +  D,  +  B, 


+  ^  +  (ro-rJ)VH 


(61b) 


46.  The  weigh.ed  residual  form  of  the  conservation  of  momentum  equations  is 
obtained  by  weighung  Equations  49  and  50  by  and  integrating  over  the  domain  Cl: 

<It  +  «  W  +  V  f  +  t'U  +  I  [?^  +  g(C  -  ar,)] 

-  H  Em  -  ^  +  riV,  ^^>a  =  0  (63) 

Applying  Gauss’s  theorem  to  the  lateral  diffusive/dispersive  terms  in  Equations  62  and 
63,  and  recalling  that  Ej^,  equals  zero  on  the  boundary,  gives  the  symmetrical  weak 

weighted  residual  form  of  the  momentum  equations: 


32 


(65) 


-  I?;  +  6(C  -  <«/)l  -  ^  +  riV,  #,>„ 

“  +  V  ^i>ii  +  ^i>n 


47.  The  GWCE  equation  is  discretized  in  time  using  a  variably  weighted, 
three-time-4evel,  implicit  scheme  for  the  linear  terms  (i.e.,  those  terms  on  the  left  side 
of  Equation  60).  Wx  and  Wy  are  treated  explicitly.  The  time  derivatives  that 
appear  on  the  right  side  of  Equation  60  are  evaluated  at  two  known  time  levels.  The 
time-discretized  GWCE  is: 


_  2  A  4.  /k-l  /k*l _ A-l 

— +  To<^  . ,  ^^i>n 

+  «i  [<g^^  .  ^>n  +  .  ^>nl 

+  «2  [<g^^  >  ^>n  +  <g^^  >  ^>nl 

+  “3  l<g^^  .  ^>n  +  <g^^  •  ^>nl 

=  .  <V^(^),  f >„ 

+  <w|,  ^>n  +  <wf, 

+  ToOr^Qn^E  +  T'otksQn**)  dP  i  =  1,  ...N 


where 

At  =  time  step 

k+1,  k,  k-1  =  future,  present,  and  past  time  levels 
ttj,  Oj,  Oj  =  time  weighting  factors 


The  time  weighting  factors  are  selected  so  that: 

OTj  +  aj  +  Oj  =  1 

o,  =  03 

Rearranging  Equation  66  gives: 
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(1  +  ^1>„ 

+  a,gAt=[<h^  ,  ^>„  +  <h^  ,  ^>(J 

+  -a^M  [<^-,  ^  <|-, 

=  2  <<‘,  h>a  +  -  1)  <C'“'.  ^i>n 

-  ajgAl2[<h^  ,  ^>a  +  <h^  ,  ^>n) 

-  .  ^>n  +  <1>|^  .  ^>nl 

+  l<r“.  ^>«  -  <t'- 

+  At<Ui‘«i'  -  ^>„  +  AKVHi"  -  ^>a 

+  At’<W!S,  ^>„  +  At=<w!;,  ^>n  -  At’F,  i 


i=l,  ...N 


where 


=  I  +  ToaiQli;!  +  ToazQS*  +  ToChQU^)  (t>i  dr  (68] 

48.  The  symmetrical  weak  weighted  residual  form  of  the  momentum  equations 
are  discretized  in  time  using  a  two-time-level  implicit  Crank-Nicolson  approximation 
for  all  terms  except  the  diffusive  terms,  which  are  treated  with  a  variably  weighted, 
two-time-level  implicit  scheme  and  the  advective,  dispersive,  and  baioclinic  terms, 
which  are  treated  explicitly: 

+  Z  <n'‘(U'‘*'  +  U''),  0i>n  -  +  V-),  *,>„ 

+  E„  [ft  <42hi-', 

,  ft  ,  ft 

=  -  z  <1:1?!^  +  -  (jS)-'’  '^On 

-  Z  4  [?!  +  sfC-  -  -  (^h)'.  *>n 


^TTk  3U'‘  .  vk  5U‘‘  2  ^  Bil  ,  ^ 

_  <U  ^  +  V  ^  ,  ^i>n  +  +  IJJ,  <f)i>^ 


i=l,  ...N 


<  — At  +  Z  +  yk),  ft>„  +<-|l'‘{u'“'  +  u^),  ^,>n 
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+  Eh!  [^1 

,  A  <^";  |(^>„ .  /*! 

=  -  J  <|^?!^‘+  ece-  -  oV“‘)l  -  (gi)-',  h>« 

-  2  <1^  1?!  +  E«‘  -  “-^ll  -  ^»n 

-  <u'‘  +  Vk  *,>„  +  <g{  +  B|,  ^,>„  i=i,  ...N  (70) 

where  /3i  and  02  are  time-weighting  factors  at  the  future  and  present  time  levels. 
These  factors  are  selected  so  that 

01  +  07  =  ^ 


Rearranging  Equations  69  and  70  gives: 

<(1  +  firjk)  Uk-J  <ln>a  -  <f'kvk‘; 

=  <(1  -  ^rjk)  Uk,  #,>„  +  <£'kvk,  ^i>„ 

-  ¥  <®  6«'‘"  -  -  (^)'‘">  '^»n 

-  At  <Uk  |2  +  ly  ’  't'i>a  ‘=k>  ■N 

<(1  +  ^rjk)  yk'l  ^i>„  +  ^  <f'kuk>i,  0i>„ 

=  <(1  -  lirjk)  yk,  4,i>„  -  ^  <t'kuk, 

-  ^!Eh!A‘  !<^‘.  |(fe>0l 

-  ¥  <|f  !?f‘+  eK"*'  -  “V"‘)l  -  (gt)'"'.  h>a 

-¥<ht^ 
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i=l,  ...N  (72) 


-  At  <U'‘  +  ]^»  ^i>n 

49.  There  are  two  differences  in  the  solution  of  Equations  67,  71,  and  72  in 
the  2DDI  option  and  the  three-dimensional  options.  In  the  2DDI  option,  the  friction 
parameter  (Cf  or  one  of  the  parameters  in  Equation  29)  is  specified  in  the  model 
input.  The  dispersive  terms  are  included  with  the  lateral  diffusive  terms  by 
eliminating  Dx  and  Dy  from  Equations  67,  71,  and  72  and  setting  E2h  =  Eah-  In  the 
three-dimensional  options,  Cf,  7,  Dx,  and  Dy  are  computed  from  the  most  recent 
internal  mode  solution.  In  flows  where  the  velocity  reverses  direction  over  the  depth, 
it  is  possible  for  the  depth-averaged  velocity  to  be  zero  while  the  bottom  stress  is 
nonzero.  In  this  case  the  drag  coefficient  computed  in  Equation  51  becomes  infinite. 
To  prevent  the  numerical  difficulties  that  this  causes,  an  upper  limit  is  set  on  the 
computed  drag  coefficient.  If  this  limit  is  exceeded,  7  and  Cf  are  set  to  zero  and  the 
bottom  stress  computed  in  the  most  recent  internal  mode  solution  is  passed  directly  to 

the  external  mode  equations.  In  the  GWCE,  —  and  ^  determined  in  the  internal 

Po  Po 

mode  solution  are  subtracted  from  W*  and  Wy,  respectively.  In  the  momentum 
equations,  and  are  subtracted  from  the  right-hand  side  of  the 
corresponding  equation.  This  modifies  the  final  terms  in  Equations  71  and  72  to 
+  HE  respectively. 

Spatial  Discretization 

50.  In  order  to  complete  the  conversion  of  the  governing  partial  differential 
equations  into  systems  of  algebraic  equations,  the  FE  method  is  applied  to  the  time- 
discretized  form  of  the  symmetrical  weak  weighted  residuwl  equations  developed  in  the 
previous  section.  Specifically,  elemental  approximations  to  the  variables  are 
substituted  into  Equations  67,  71,  and  72,  the  elemental  equations  are  summed  over 
the  global  domain,  and  the  required  degree  of  inter-element  functional  continuity  is 
enforced.  Interpolating  basis  with  at  least  C®  functional  continuity  are  required  to 
discretize  most  of  the  dependent  variables.  Departures  from  this  are  noted  below. 

51.  In  all  linear  terms,  surface  elevation,  velocities,  and  depth  are 
approximated  over  each  element  as: 
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(73a) 


<‘‘5  E  <t'i 

j  *  1  ■'  •* 


^  E  Uf 

j  *  i 


V'‘  5  E 
j  =1 

J1  e  1 


(73b) 
(73c) 

h  2!  E  h:  (73d) 

j  .  1  J  * 

where  equals  the  number  of  nodes  per  element.  In  nonlinear  terms  and  certain 

forcing  terms,  the  entire  term  may  be  interpolated  over  the  element  as  described 
below. 

52.  The  nonlinear  and  forcing  terms  in  the  GWCE  are  approximated  as 
follows. 

a.  The  Coriolis  parameter  and  the  fluxes  in  the  Coriolis  term  are 
appro.'dmated  by: 

f'’‘(UH)>‘  ~  ”e^  (f'UH)f  4>-,  (74a) 

i  =1  *  * 


f/k(vH)k  !v  (f'VH)^  4>. 
j  « 1 

fe.  The  finite  amplitude  component  of  the  free  surface  gradient  is 
approximated  by: 

«7)k  ?  “e  (OJ 

j  *  1 


(74b) 


(75) 


c  The  combined  barometric  pressure  and  Newtonian  tidal  potential  term 
is  approximated  by: 

(5?i  - 


Vog  "  -  j  M  Vog 

The  total  depth  factor  in  this  term  is  evaluated  using  an  L: 
approximation: 

2;  ^  “e^  H5 

-  nei  j  J 

d-  The  surface  stress  terms  are  approximated  b3'; 


Ue  1 

(Isi)k  ^  6. 

Vo  -  j  =1  Vo 


^77) 

(78a) 

(78b) 
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e.  The  bottom  stress  and  Tq  terms  are  approximated  by; 

(rik-  !  “e’  [(rj  -  ro)HUlf  A 

j  *  1  ^ 


(79a) 


(Ti'‘-  To)H'‘V'‘  5  E  l(r;  -  r„)HV)}  (79b) 

j  *  I 

In  the  three-dimensional  model  options,  if  the  computed  friction 
coefficient  exceeds  the  maximum  allowable  value,  the  bottom  stress 
terms  are  approximated  directly  by: 


k  Del 

Iki  ~  JJ  {'Ibs)k  A. 
Po  =  jTi  Vo 

N  *^1;^  (Ibi)k  A 

Po  =  jTi  Vo  h 


(80a) 


(80b) 


f.  The  dispersive  terms  are  broken  up  into  their  components  Duu,  Duv 
and  Dw,  (defined  in  Part  II),  and  discretized  as; 


S  “e_  ^ 

I^Dvv)'-  S  ^ 


(81a) 


(81b) 


(81c) 


(81d) 


g.  The  baroclinic  terms  are  not  included  in  either  ADCIRC-2DDI  or 
ADCIRC-3DL.  Therefore  the  discretization  of  these  terms  is  not 
considered  here. 

h.  The  velocities  that  multiply  the  time  derivative  components  of  the 
non-conservative  advective  terms  are  approximated  using  L2 
interpolation: 


^  =  - 


1  “el  , 
el  j  =  t  •' 


(82a) 


v‘  =  E  5^  ”e  VJ  (82b) 

The  free  surface  elevations  that  appear  in  these  terms  are 
approximated  using  the  standard  Co  approximation.  Equation  73a. 
The  spatially  differentiated  components  of  the  non-conservative 
advective  terms  are  approximated  by 
mi  ,  ,  ®el  . 


(UH  §)‘  =  (UH)J,  ,E  U}  ^ 


(83a) 
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(VH  1^)“  !  (VH)5,  U}  ^  (83b) 

(UH  5  (UH)5i  ”s‘  V8  ^  (83c) 

(VH  |^)‘  S  (VH)J,  vy  ^  (838) 

where 

f  ^  6  1 

(UH)Si  H  ^  E  (UH)}  (84a) 

“el  j  =  I 
1  Q  e  1 

(VH)5,  e  i  .E  (VH)8  (84b) 


53.  Substituting  the  approximations  in  Equations  73  -  84  into  Equation  67 
and  summing  over  the  elements  gives: 

He  1 


H 

E 

e  1  »l 


+ 

+ 


j?.  [d  + 

<.igAl2[<  ||i,  ^>0,,+  <^S  ^>nj 

^  ^>0.,  +  «i’‘ 

-  “i8^d[<^E^  h„0„<y  ^>n„+  j  M.<i 

-  h.0.<S-  <”s_ 

^  «4-  ^>0,  +  «S-  f.  ^>«J 

At  l<US.(fl  -  <f-)  4>j,  +  <V‘,«5  -  <?-')  fi>„J 

-  Att[<(UH)5.u5  gi>„^,+  <(VH)S.US 

+  <(UH)S.VS[  +  <(VH)‘,V} 

+  Attl<(f'VH)8  -  <(rUH)f  ^>„J 

-  ^  1<(05  ^>8.,  +  <«’)?  f  >»,.i 

-  «AttHj.l<(^  -  a,))  ^  >„,,+  <(^  -  m,))  |i.  fi>„.,l 


+ 

+ 
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+ 

-  At'{<l(rj  -  r.)HU|5  +  <[(rj  -  r.)HV];  0^, 


where 

M  =  the  total  number  of  elements 

flgj  =  the  elemental  domain  for  the  element,  el. 


Equation  85  can  be  rewritten  as; 


M 

E 

e  1  *t 


nei 

E 

j  *1 


[(1  + 


=  "e‘  [2m(|)<J  +  (2^  -  -  a,gAt^M(?)<J 

+  IM^SlCCf  -  Cf")  +  -  <}-)j 

-  At>(M(t>(UH)‘.uf  +  M(f(VH)5,u}  +  m!?(uh)5,v}  +  m(?)(vh)5,v;i 


+  At*(M(p(f'VH)J  -  M(p(f'UH)Jl  -  MjP(Oj 


-  At*{M(T)[(r;  -  r.jHU]}  +  M(f)[(rJ  -  rojHV]}) 

-  At>[M(pD„„J  +  M(f)D„vf  +  M(pD„,J  +  M(pDv,{]] 


-  At’F, 
i=l,  ...N 


where 


Del 

+  <  E 

m  *  1 


Mi?)  ^  KK  t- 

Mi?)  s 


(86) 

(87a) 

(87b) 

(87c) 
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mW  .  <^.  ^>0^, 

mIP  ^  <^i>  ^>0., 

mSP  =  <^j. 

Note  that  the  elemental  matrices,  Mi})  ,  mW,  m}})  .  Mif,)  and  Mif)  are 
symmetrical  and  that  M^p,  and  M^p  are  non-symmetrical. 

54.  The  fully  discretized  GWCE  can  be  written  in  a  compact  form  as: 

E  {Cf  1}  =  E  i=l,  ...N 

c  1  “1  j  « 1  ■*  e  1  »1 


(87d) 

(87e) 


(87g) 

(87h) 


where 


Mff®®  s  (1  +  ^)Mi})  +  a^At>Mif)  +  56^  m}})  ( 

pewcE  ^  ”1;'^  [sMi})(}  +  -  l)Mi})<}-'  -  a,gAt’Mif)<} 

-  arfAt=Mif)<}-‘  +  Mif)<}-‘ 

+  At  (Mi})uk,(C}  -  <}■•)  +  Mif)vk,«}  -  <}-)) 

-  At^Mif)(UH)‘,U;  +  Mif)(VH)iiU}  +  Mif)(UH)5,V5  +  Mif)(VH)5iV5l 
+  At*IMi})(f'VH);  -  Mif)(f'UH);]  -  Mif)(OJ 

-  8A‘’Mif)(^  -  +  At=[MiP(^)}  +  Mif)(I“)}l 

-  AtHMi})[(n  -  t„)HU)}  +  mWKtJ  -  r,)HVl}} 

-  At’(Mi})D„„5  +  Mif)D„v}  +  Mif)D„,5  +  Mif)D„}l]  -  At>Fi 

i=l,  ...N  ( 

In  the  three-dimensional  options,  if  the  computed  friction  coefficient  exceeds  the 

CWCE 

maudmum  allowable  value,  the  friction  term  in  the  right  side  load  vector  Pj  is 
slightly  different: 
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pGWCE  ^  “g‘  +  (2^  -  -  o^At“M(?)c5 

+  At  [MipUk,((:V  -  <lf->)  +  M(®)v5,«j[  -  C|f')l 

-  At’[M(t)(UH)‘,Uj[  +  Mlf)(VH)‘,U}  +  M®(UH)‘,V5  +  M(^)(VH)k,V}l 

+  At’IM^JWjf  -  Ml?){roH)J)  -  mSJ) (05 

-  gAtm)(^  -  <»,)5Hk,  +  At=lM(J)(lM)k  +  m(?)(2«)}) 

-  At“{M(J)[(^)5  -  r,(UH)5]  +  m(®)[(J^)5  -  t„(VH)51} 

-  At>(Mjt'D„„}  +  m(®)d„,5  +  mWd„,5  +  m(^)Dv,5]]  -  At»Fi 

i=l,  ...N  (91) 

55.  Global  assembly  and  enforcement  of  the  functional  continuity 
requirement  leads  to  the  following  global  system  of  equations: 

E  [  ]  {  8Ci*‘  }  =  {  }  j=l.  •••N  (92) 

j  M  ■*  J 

where 

g^GWCE  _  global  banded  system  matrix 

gpGWCE  _  global  load  vector 

=  the  global  nodal  elevation  vector 

56.  The  fully  discrete  form  of  the  momentum  equations  is  obtained  from  the 
time-discretized  symmetrical,  weak  weighted  residual  form  of  the  momentum  equations, 
Equations  71  and  72,  as  follows. 

g,.  The  local  acceleration  terms  are  interpolated  using  Equations  73b  and 
73c. 


The  friction  terms  are  approximated  by: 

j  =  i 

(93a) 

j  *  1  •' 

(93b) 
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(93c) 


Del 

j  •  i  ■*  ■' 


(93d) 


In  the  three-dimensional  model  options,  if  the  computed  firiction 
a)efGcient  exceeds  the  maximum  allowable  value,  the  bottom  stress 
terms  introduced  on  the  right-hand  side  of  the  momentum  equations 
are  approximated  by: 


ttel 


(94a) 


(94b) 


c.  The  Coriolis  terms  are  approximated  by: 
^  €  1 

f.kukM  ^  2  f/kuk*i  ^ 

j  *  1  '*  ■'  •* 

nel 

f/kvk»i  ~ 

j  =  1  ^  ^  ^ 

nei 

f/kyk  2;  £  f/kuk  ^ 

j  =  1  •' 

nei 

f/kyk  N  S  f/kyk  ^ 
i  «i  ■' 


(95a) 

(95b) 

(95c) 

(95d) 


d.  The  lateral  diffusive/dispersive  terms  are  approximated  using 
Equations  73b  and  73c  for  velocity  and  Equation  77  for  total  depth. 

e.  The  barometric  pressure  and  Newtonian  tidal  potential  are 
interpolated  using  Equation  76. 

f.  The  surface  elevation  is  approximated  using  Equation  73a. 

g.  The  surface  stresses  are  evaluated  as: 


n  e  1 

/'^SX'ik*!  ni  V 

W 

(96a) 

n  e  1 

^  J?, 

(96b) 

n  e  1 

(96c) 

^  €  2 

(96d) 

h.  The  advective  terms  are  approximated  by: 
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r  =  V,V£;  U5  ^  (97b) 

u^  t  USi^E  Vf  ^  (97C) 

V"'  W  ^  VS, “2,  Vj  ^  (97d) 

where  U^i  and  V^j  are  defined  in  Equation  82. 

i.  The  dispersive  terms  are  broken  up  into  their  components  Duu,  Duvi 

and  Dw,  (defined  in  Part  II),  and  discretized  as: 

(h  f  “"“f  ^ 

(h  f ”  (hIi)''”?.  “-i  ^  (»*'>) 

(h  f®*') 

(h  f s  (0^:)““?, 

where  Hei  is  defined  in  Equation  77. 

j.  The  baroclinic  terms  are  not  included  in  either  ADCIRC-2DDI  or 


ADCIRC-3DL.  Therefore  the  discretization  of  these  terms  is  not 
considered  here. 

57.  Substituting  the  approximations  in  Equations  93  -  98  into  Equations  71 
and  72  and  summing  over  the  elements  gives  the  discrete  system  of  equations: 
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-  A‘l<U^iUj  #i>n.,  +  <V5,U5  *,>„J 

-  Atl<(jL)kD„„k  +  <(jL)kD„,k  01>()j] 

i=l,  ...N 


J..  %  ^  <f'jUi*Vi. 

=  <(1  -  t^rij)  Vj0j. 

-  ^jEhiAt  [<VkH|i  |^^)>(i^,+  <VjH|,  ^^)>nj 

-  +  Ci“  -  “^“)  +  (?!i  +  <i  -  ^)1  *i>»., 

+  l^t(^)j“  +  (^)il  *l>nel 

-  At(<U5iVj[  4k,>„^,+  <V‘,Vk  ^j>„J 

-  At[<(g^)'‘D„vJ  «i>n^j  +  <(h^)'‘Dvv5  ^i>nj] 


i=l,  ...N 


Equations  99  and  100  can  be  rewritten  as: 

“  ®el  r  At  . .  m  .  .  At  .  I 


_E__  E  [(1  +  ^4j)M(!)uf  •■  -  +  ^,Et,AtMS?)u5' 

=  (1  -  firjJjM^pUV  +  Aif.kjjWvk  -  ,3,Et,AtM!3)u5 

-  +  (5“  -  a^-)  +  (^  +  C$  -  o»^)l 

+  +  vJimWuJ) 

-  At(^)'‘lM(pD„„f  +  mWDuv}]]  i=l.  ...N 


(100) 


(101) 


I?.  7^^il)MipV5-  +  fif'jMSjJuJ*'  +  ^.E,,AtM(fvJ- 
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58.  The  M(p  matrices  on  the  left  side  of  Equations  101  and  102  and  in  the 

first  two  terms  on  the  right  side  of  these  equations  are  lumped  so  that  all  elements 
are  added  onto  the  diagonal.  The  matrices  on  the  left  side  of  Equations  101 

and  102  are  decomposed  into  diagonal  and  non-diagonal  matrices.  The  non-diagonal 
portion  of  is  moved  to  the  right  side  of  the  equations.  These  operations  give: 


-  +  <j“  - 

+  -  At(US,M(I)u{  +  Vj,MSf)up 

-  At(Hl)'‘lMWD„„}  +  MWD„,y]]  i=l,  ...N  (103) 

and 

E  [(1  + 

C  1  *1  J  ”  ^  ^ 

=  (1  -  |irt5)Ml5‘')vy  -  -  E^5At(^,M(f'>)vV->+  ^jMS^V}) 

-  +  (}•'  -  +  (^  +  Cf  -  ao})] 

+  +  (^)jl  -  +  v5,m(«)v}) 

-  Al(3L)k[M5pD„vf  +  Mj?)Dvv5l]  i=l,  ...N  (104) 

where 

=  the  diagonally  lumped  elemental  matrix 

-  the  diagonal  portion  of  the  elemental  matrix 
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Mlf*)  =  the  non-diagonal  portion  of  the  elemental  matrix 


59.  The  fully  discretized  momentum  equations  can  be  written  in  compact  form 


M  Uel 


E  S  [Mtn  -  [Mfi  ]  {VVn  =  E  {Pt  } 

1=1  j  =1  J  ‘  J  ■*  J  J  ->  el=l 


i=l,  ...N 


M  Uel 


E  E  [Mfi  1  ]  {VV‘}  =  2  {P;  }  i=l,  ...N 

1  *1  j  =1  •'  ■*  ■'  •*  e  1  *1 


where 


Mjf  =  (1  +  +  /3.Et,AtM(f ) 


j^2me 

ij  2  J 


pI»E  ^  "£‘  ^ 

-  EuAt(/?,M(f  +  ^kHWuf) 

-  8|iM(T)[(^'  +  CV“  -  m/J-)  +  (^  +  Cf  -  0^)) 

+  +  (gt)S  -  ^‘(''^1  +  V‘,M}?)uf) 

-  At(gL)'‘(M(I)D„„J  +  m(^)d„,J]]  i=l,  ...N 

pT"'=  =  ”e‘  ^(1  -  |iTi5)M(‘‘'V5  - 

-  EhjAt(/J,M(f  "Vf  + 

-  +  c;-‘  -  0^")  +  (^  +  C{  -  0-^)1 

+  +  (^)j!  -  luSiwSPvJ  +  v‘,m(?)v5] 

-  At(jjL)k[M(pD„v}  +  m(?)d,vSi]  i=l,  ...N 


j. 

+  -r* 


i=l,  ...N 


(105) 


(106) 


(107a) 

(107b) 


(107c) 


(107d) 


In  the  three-dimensional  model  options,  if  the  computed  friction  coefficient  exceeds  the 
maximum  allowable  value,  the  bottom  stress  terms  appear  explicitly  on  the  right  side 

XME  YMF 

of  the  momentum  equations  and  therefore  are  included  in  Pj  and  P^  : 


j7j  ^  u 
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-  ■  “’^>1 
+  +  (^)5]  -  At(Uki  M(pU5  +  Vk,M(?)uy) 

-  At(5^)'‘lM(I)D„„5  +  m(?)d„v}|  - 

p^E  ^  “|,>  -  ^m(}‘')u} 

-  Et,At(/3,M(f”)vf>  + 

-  +  <l"  -  “^*')  +(?!+<?- 

+  +  (^)5|  -  At(Uj,M(J)v5  +  V‘,M(pv5) 

-  ^‘(h^)‘1m1?D»vJ  +  M(f)Dv.;]  -  AtM(?)(^);' 


(108a) 


(108b) 


60.  Global  assembly  and  enforcement  of  the  C°  functional  continuity 
requirement  gives  the  following  systems  of  equations: 

N 

S  [eM|f]  {gUf ‘}  -  {gM?f]  {gV^*‘}  =  i=l,  ...N 

j  j  j 

N 

E  {gUj*‘}  +  [gM{^®]  {gyf =  {gpf**^}  i=l,  ...N 

where 

gM^^®,  =  global  diagonal  system  matrices 

gpXME^  gp^E  _  giQij^  right-hand-side  load  vectors, 

gUj*^  gVj*‘  =  global  velocity  vectors  in  the  x  and  y  directions 


(109) 


(110) 


Solution  StrateE 


61.  The  horizontal  discretization  for  ADCIRC  has  been  implemented  with 
three-node  linear  triangles  and  four-node  bilinear  quadrilaterals.  The  triangle  element 
provides  a  maximum  degree  of  flexibility  and  is  extremely  cost-effective  on  a  per-node 
basis  for  long  wave  computations.  All  elemental  matrices  through  are 

integrated  using  a  numerical  quadrature  rule  that  is  specified  with  the  input  data.  A 
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four-point  Gaussian  quadrature  rule  integrates  the  elemental  matrices,  Equation  87, 
exactly  (Connor  and  Brebbia  1976).  However,  for  most  applications,  a  three-point 
Gaussian  quadrature  rule  appears  to  be  sufQcient.  The  elemental  are 

computed  once  and  then  stored  for  use  during  the  time-stepping  operations. 

62.  The  GWCE  is  solved  first.  The  global  system  matrix  for  the  GWCE, 

SMjj  ,  is  time-independent  and  is  therefore  assembled  and  LU  decomposed  only 

once.  ®Mij  has  a  banded  structure  with  a  band  width  that  is  dependent  on  the 

node  numbering  of  the  grid.  Prior  to  running  ADCIRC,  the  node  numbering  should 
be  optimized  to  minimize  the  maximum  difference  in  the  global  node  numbers 

GWCE 

associated  with  each  element  in  the  grid.  The  right  side  load  vector,  sp.  ^  ig 

efficiently  updated  every  time  step  since  all  elemental  matrices  have  been 
pre-computed.  Flux-spedfied  boundary  conditions  have  been  incorporated  into  the 
load  vector  by  the  model  formulation  and  therefore  require  no  additional  equation 
manipulation.  Elevation-specified  boundary  conditions  are  incorporated  into  the 
system  matrix,  ,  by  zeroing  out  rows  corresponding  to  boundary  nodes  with 

specified  elevations  and  placing  a  value  of  unity  onto  the  diagonal.  The  boundary 
condition  values  are  then  stored  into  the  appropriate  location  in  The 

equations  associated  with  elevation-specified  boundary  conditions  are  multiplied  through 
by  a  constant  to  ensure  that  the  modified  global  matrix  has  a  good  condition  number. 
The  modified  system  of  equations  (i.e.,  Equation  92  modified  to  include  the  elevation- 
specified  boundary  conditions),  is  then  solved  for  elevation  at  the  new  time  level,  k+1. 

63.  The  momentum  equations  are  solved  second  and  use  the  elevation  values 

at  time  level  k+1  computed  from  the  GWCE.  The  global  system  matrices  for  the 
momentum  equations,  and  ),  are  time-dependent  and  therefore  need  to 

be  reevaluated  at  every  time  step.  However,  since  these  matrices  are  diagonal,  matrix 
evaluation  and  decomposition  are  very  economical.  Specified  normal  flux  boundary 
conditions  are  incorporated  into  Equations  109  and  110  by  reorienting  the  x  and  y 
equation  pairs  that  correspond  to  the  specified  flux  boundary  nodes  into  a  locally  (for 
each  node)  normal/tangential  coordinate  system.  The  i3oriented  equations  are  then 
replaced  by  the  corresponding  specified  normal  flow  boundary  condition  values  (Wang 
and  Connor  1975;  Gray  1984). 

64.  The  right  sides  of  Equations  109  and  110  are  dependent  on  Uj, 

and  Vj,  which  are  all  known  quantities,  and  on  and  (because  of  the  lateral 
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closure  model).  Therefore,  these  equations  must  be  solved  iteratively  for  velodti^  at 
the  new  time  level,  k+1.  and  sp.  are  updated  each  iteration  with  the  new 

values  of  Uj**  and  Vj*‘  until  a  specified  convergence  criteria  has  been  reached.  When 

is  zero,  no  iteration  takes  place. 

65.  When  a  three-dimensional  option  is  used,  the  external  mode  solution 

depends  on  the  internal  mode  solution  through  Duu)  Duv>  Dw,  Cf  and  7  [or  if  the 
drag  coefficients  exceed  the  maximum  allowable  values  on  and  These 

quantities  are  computed  at  each  internal  mode  time  step  and  assumed  to  be  constant 
in  time  for  subsequent  external  mode  time  steps.  If  the  external  mode  solution  and 
the  internal  mode  solution  are  evaluated  at  the  same  model  time,  the  external  mode 
solution  is  evaluated  first.  The  updated  surface  elevations  and  depth-averaged 
velocities  are  then  used  in  the  internal  mode  solution.  This  solution  sequence  requires 
the  specification  of  initial  values  i  c  Duu>  Duv,  Dw,  Cf,  and  7  as  input  parameters  for 
the  external  mode  solution. 

Fourier  Properties  of  the  External  Mode  Solution 

66.  Fourier  analysis  characterizes  the  damping  and  phase  propagation  properties 
of  a  numerical  solution  in  relation  to  the  corresponding  analytical  solution.  Although 
it  is  typically  applied  to  the  one-dimensional  form  of  the  shallow-water  equations  and 
a  constant  bathymetric  depth  is  usually  assumed,  the  results  give  a  good  indication  of 
how  a  circulation  model  will  behave  in  a  more  general  two-dimensional,  nonlinear  field 
application.  They  also  allow  inter-comparisons  with  other  discretization  strategies. 
Procedures  for  applying  Fourier  analysis  to  the  shallow-water  equations  are  described 
by  Pinder  and  Gray  (1977)  and  Lynch  (1978). 

67.  The  discrete  form  of  the  ADCIRC  2DDI  governing  equations  has  been 
Fourier  analyzed.  These  results  are  presented  below  along  with  results  from  the 
Fourier  analyses  of  several  other  numerical  solution  schemes  for  the  shallow-water 
equations.  All  of  the  other  numerical  schemes  that  were  considered  use  primitive 
formulations  of  the  shallow-water  equations  (as  opposed  to  the  generalized  wave- 
continuity  formulation  used  in  ADCIRC).  The  schemes  include  a  finite  element 
solution  using  linear  elements  (PEFE)  (Wang  and  Connor  1975;  Westerink,  Connor 
and  Stolzenbach  1987,  1988),  a  second  order,  non-staggered,  finite  difference  solution 
(PENSFD),  and  a  second  order,  staggered,  finite  difference  solution  (PESFD)  (Hansen 
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1956;  Leendertse  1967).  A  second  order,  Crank-Nicholson  scheme  was  used  to 
integrate  the  PEFE,  PENSFD,  and  PESFD  in  time.  As  described  previously, 

ADCIRC  uses  a  three-time-level  scheme  for  the  GWCE  (oi  =  0.35,  oj  =  0.30  and 
as  =  0.35),  and  a  Crank-Nicholson  scheme  for  the  momentum  equations.  The  bottom 
friction  coefficient  (Equation  46)  in  each  model  was  specified  as 

r*  =  0.8  If  /g5~/A  (111) 

where  X  is  the  wavelength  of  the  Fourier  component. 

68.  The  modulus  of  the  propagation  factor  indicates  the  ratio  of  the  numerical 
amplitude  to  the  analytical  amplitude  during  the  propagation  of  one  wavelength.  The 
phase  of  the  propagation  factor  indicates  the  phase  lag  or  lead  a  given  wavelength 
experiences  during  one  period.  Figure  1  pr^ents  the  modulus  and  phase  of  the 
propagation  factor  for  the  PEFE  and  PENSFD  schemes.  Comparisons  are  shown  for 
Cr  =  0.1,  0.5,  1.0,  and  2.0  where  Cr  is  the  Courant  number  based  on  wave  celerity. 

At  is  the  time  step,  and  Ax  is  the  gnd  spacing.  For  increasing  Cr,  both  the  PEFE 
and  the  PENSFD  solutions  have  less  damping  than  the  analytical  solution  for  low 
ratios  of  A/ Ax.  Neither  solution,  regardless  of  Cr,  propagates  energy  at  the  shortest 
resolvable  wavelength,  A  =  2Ax.  This  characteristic  of  PEFE  and  PENSFD  solutions 
accounts  for  the  severe  2Ax  numerical  noise  problems  encountered  using  these  schemes. 

69.  Figure  2  presents  the  modulus  and  phase  of  the  propagation  factor  for  the 
PESFD  scheme  and  the  generalized  wave-continuity  equation  finite  element 
(GWCEFE)  scheme  used  in  ALCIRC.  For  low  ratios  of  A/Ax,  both  schemes  provide 
less  damping  than  the  analytical  solution  and  show  poorer  phase  propagation  behavior 
as  Cr  increases.  For  a  fixed  Cr  and  A/ Ax,  the  PESFD  scheme  has  slightly  better 
damping  characteristics,  while  the  GWCEFE  scheme  has  better  phase  propagation 
characteristics.  At  low  Cr,  the  GWCEFE  solution  leads  the  analytical  solution.  As  Cr 
increases,  the  GWCEFE  phase  propagation  factor  swings  through  a  zero  value 
(corresponding  to  perfect  phase  behavior)  and  then  develops  a  phase  lag.  This 
indicates  that  there  will  be  a  local  minimum  in  the  time  convergence  curve  with 
optimal  accuracy  being  achieved  at  Cr  »  0.5. 

70.  The  primary  difference  between  numerical  solutions  using  PEFE  and 
PENSFD  schemes  and  numerical  solutions  Uoing  GWCEFE  and  PESFD  schemes  is 
that  the  latter  schemes  propagate  energy  at  A  =  2Ax.  Propagation  of  2Ax  waves 
corresponds  to  a  non-folded  dispersion  relationship  and  prevents  two  responses  from 
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MODULUS  OF  PROPAGATION  FACTOR  MODULUS  OF  PROPAGATION  FACTOR 


a.  Cr  =  0.1 


J/Ax 


b.  Cr  =  0.5 


Figure  1.  Modulus  and  phase  of  the  propagation  factor  for  PEFE  and  PENSFD 
solutions  (Continued) 
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MODULUS  OF  PROPAGATION  FACTOR  MODULUS  OF  PROPAGATION  FACTOR 


VAX  VAX 

a.  Cr  =  0.1 


^Ax  VAX 


b.  Cr  =  0.5 


Figure  2.  Modulus  and  phase  of  the  propagation  factor  for  GWCEFE  and  PESFD 
solutions  (Continued) 
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MODULUS  OF  PROPAGATION  FACTOR  MODULUS  OF  PROPAGATION  FACTOR 


developing  to  a  single  forcing  frequency,  i.e.,  one  physical  response  at  the  forcing 
wavelength  and  one  numerical  response  at  a  wavelength  near  2Ax  (Platzman  1981). 

As  a  consequence,  GWCEFE  and  PESFD  schemes  do  not  have  the  severe  2Ax  noise 
problem  of  the  PEFE  and  PENSFD  schemes. 

Convergence  Properties  of  the  External  Mode  Solution 

71.  In  order  to  verify  the  accuracy  of  the  external  mode  solution  of  ADCIRC 
and  to  establish  convergence  properties  in  space  and  time,  ADCIRC-2DDI  was  applied 
to  a  modified  form  of  the  quarter  annular  test  problems  originally  developed  and 
appMed  by  Lynch  and  Gray  (1978,  1979)  and  Gray  and  Lynch  (1979).  These  two- 
dimensional,  variable-depth  test  problems  were  developed  to  give  insight  into  a 
numerical  scheme’s  2Ax  oscillations  and  its  ability  to  propagate  longer  physical  waves. 
The  original  geometry  and  bathymetry  of  Lynch  and  Gray  (1978,  1979)  were  modified 
as  follows.  The  arc  of  the  annulus  was  increased  to  135  deg*;  the  inner  radius  was 
decreased  to  125,000  ft;  the  outer  radius  was  increased  to  650,000  ft.  The  resulting 
geometry,  with  three  land  boundaries  and  one  open  ocean  boundary,  is  shown  in 
Figure  3.  A  linearly  varying  bathymetry  was  used  that  increased  from  50  ft  at  the 
inner  radius  to  260  ft  at  the  outer  radius  and  a  quadratically  varying  bathymetry  was 
used  that  increased  from  50  ft  at  the  inner  radius  to  1,352  ft  at  the  outer  radius. 
These  modifications  accomplish  two  things.  First,  the  modified  domains  are  more 
representative  of  a  coastal  region  that  extends  to  near  or  beyond  the  Continental  Shelf 
break.  (In  fact,  the  geometry  and  bathymetry  are  idealized  approximations  to  the 
New  York  Bight.)  Second,  the  numerical  difficulty  of  the  test  problems  is  increased. 

72.  A  sequence  of  four  discretizations  was  considered;  a  6-  by  8-node 
discretization  (Ar  =  105,000  ft),  an  11-  by  15-node  discretization  (Ar  =  52,500  ft),  a 
21-  by  29-node  discretization  (Ar  =  26,250  ft),  and  a  41-  by  57-node  discretization 
(Ar  =  13,125  ft).  These  are  shown  in  Figure  4.  Grids  consisting  of  linear  triangles 
and  of  bilinear  quadrilaterals  were  tested  and  gave  very  similar  results.  Only  the 
bilinear  quadrilateral  results  are  presented  here.  For  each  grid,  five  different  time 
steps  were  applied:  At  =  T112/8,  At  =  TM2/I6,  At  =  TM2/32,  At  =  TM2/64,  and 

At  =  Tii2/128  where  T11I2  is  the  M2  tidal  forcing  period  equal  to  44,712  seconds. 
ADCIRC-2DDI  was  run  in  its  linear  mode  with  an  M2  forcing  frequency.  Therefore, 


*A  table  of  factors  for  converting  non-SI  units  of  measurement  to  SI  (metric)  units  is 
presented  on  page  6. 
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Figure  3 


650.000  ft. 


Test  problem  geometry 
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a.  6-  by  8-node  grid 


b.  11-  by  15-node  grid 


Figure  4.  Grids  used  for  the  test  problem  (Continued) 
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d.  41-  by  57-node  grid 
Figure  4.  (Concluded) 
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the  theoretical  model  response  should  have  included  only  an  Mj  wave.  The  resolution 
of  the  Mj  wave  provided  by  the  sequence  of  grids  varied  between  17  and  312  nodes 
per  wavelength  for  the  linearly  varying  bathymetries  and  between  17  and  703  nodes 
per  wavelength  for  the  quadratically  varying  cases.  Thus  the  M3  wave  was  always 
well-resolved.  €r  varied  between  0.13  and  39  for  the  linearly  varying  bathymetries 
and  between  0.13  and  88  for  the  quadratically  varying  cases.  ADCIRC-2DDI  is 
unconditionally  stable  in  its  linear  mode  and  therefore  permits  the  use  of  Cr  >  1. 

73.  All  cases  were  forced  at  the  open  ocean  boundary  using  ^  =  1.0  8in(u;|^  t) 

where  0;^^=  2t/TiI3  is  the  Ms  forcing  frequency.  All  other  forcing  mechanisms  (i.e., 

tidal  potential,  free  surface  wind  stress  and  atmospheric  pressure  gradients)  were  set  to 
zero.  The  Coriolis  and  advective  terms  were  also  neglected.  The  bottom  friction 
coefficient  was  set  to  r*  =  0.0001  and  the  value  of  Tq  =  0.0001.  All  total  depths 


were  set  equal  to  the  depth  to  the  geoid. 

74.  The  computations  were  hot-started  using  the  analytical  solution  for  the 
specified  geometry,  bathymetry,  and  friction  coefficient.  The  computations  were  then 
run  for  10  tidal  cycles  to  allow  a  dynamically  steady-state  numerical  solution  to 
develop.  The  elevation  and  radial  velocity  solutions  at  each  node  were  recorded 
during  the  eleventh  tidal  cycle  and  were  Fourier  decomposed.  Typical  results  are 
shown  in  Figure  5  for  the  sequence  of  runs  using  the  coarsest  grid  and  the  linearly 
varying  bathymetry.  The  figures  compare  the  exact  analytical  solution  to  the 
maximum  and  minimum  ADCIRC-2DDI  solution  for  all  nodes  at  the  same  radius. 
These  plots  indicate  that  there  are  no  spurious  2Ax  modes  in  either  the  radial  or 
angular  directions. 

75.  Error  measures  were  calculated  from  comparisons  between  the  harmonically 
decomposed  numerical  solutions  and  the  analytical  solutions.  These  were  defined  as: 


(113a) 

(113b) 

(113c) 

(113d) 
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R  Xio^ft)  R  XlO^ft) 

a.  At  =  TM2/128 

Figure  5.  Elevation  and  radial  velocity  components  on  the  6-  by  8-  node  grid  with  linear  bath3rmetry 
(Sheet  1  of  5) 
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FigTJie  5.  (Sheet  5  of  5) 


where 


Np  =  the  number  of  nodes  within  the  grid 

=  amplitude  of  the  sine  component  of  the  analytical  elevation  solution  at  node  i 

bf  =  amplitude  of  the  cosine  component  of  the  analytical  elevation  solution  at  node  i 


n 

a^  =:  amplitude  of  the  sine  component  of  the  numerical  elevation  solution  at  node  i 


amplitude  of  the  sine  component  of  the  analytical  radial  velocity  solution  at 
node  i 

amplitude  of  the  cosine  component  of  the  analytical  radial  velocity  solution  at 
node  i 

amplitude  of  the  sine  component  of  the  numerical  radial  velocity  solution  at 
node  i 


amplitude  of  the  cosine  component  of  the  numerical  radial  velocity  solution  at 
node  i 


These  error  measures  represent  the  absolute  errors  in  the  sine  component  of  the 
elevation  solution  (El),  in  the  cosine  component  of  the  elevation  solution  (E2),  in  the 
sine  component  of  the  radial  velocity  solution  (E3),  and  in  the  cosine  component  of 
the  radial  velocity  solution  (E4). 

76.  A  summary  of  the  error  measures  computed  for  all  of  the  test  runs  is 
presented  in  Table  2.  The  error  measures  are  plotted  against  Cr  (the  average  value 
for  a  given  grid)  for  the  linear  bathymetry  test  cases  in  Figures  6  and  7  (Figure  7  is 
a  blow-up  of  the  low  Cr  range  in  Figure  6),  and  for  the  quadratic  bathymetry  test 
cases  in  Figures  8  and  9  (Figure  9  is  a  blow-up  of  the  low  Cr  range  in  Figure  8). 

All  errors  show  good  spatial  convergence;  i.e.,  the  more  refined  the  grid,  the  lower  the 
error  at  any  Cr-  In  time,  the  errors  decrease  as  Cr  decreases,  until  Cr  =  0.9  -  1.75 
for  the  linear  bathymetries  and  Cr  =  3.5-7  for  the  quadratic  bathymetries.  A  well- 
defined  local  error  minimum  exists  for  all  grids  within  these  Courant  ranges  for  both 
the  sine  and  cosine  components  of  the  elevation  and  radial  velocity  solutions.  This 
local  error  minimum  occurs  because  the  phase  of  the  propagation  factor  changes  from 
a  phase  lead  to  a  phase  lag,  passing  through  a  region  of  almost  perfect  phase 
behavior,  near  Cr  »  0.5  (see  Figure  2  and  associated  discussion).  Figures  6-9 
suggest  that  the  optimal  behavior  occurs  at  somewhat  higher  values  of  Cr.  These 
figures  were  plotted  using  the  average  value  of  Cr  for  a  given  grid.  However,  the 
primary  errors  are  generated  in  the  shallow  portions  of  the  domain.  If  the  Cr  is 
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Figure  6.  Convergence  properties  of  the  linear  bathymetry  test  cases 
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Figure  7.  Detail  of  the  convergence  properties  of  the  linear  bathymetry  test  cases  at  low  C; 
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Figure  8.  Convergence  properties  of  the  quadratic  bathymetry  test 
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convergence  properties  of  the  quadratic  bathymetry  test  cases  at  low  C 


adjusted  to  account  for  this,  the  optimal  range  of  values  changes  to  Cr  =  0.52  -  1.07 
for  the  linearly  varying  bathymetries  and  Cr  =  1-07  -  2.14  for  the  quadratically 
varying  bathymetries. 

77.  It  is  concluded  that  the  external  mode  solution  used  in  ADCIRC  has 
excellent  numerical  properties.  There  are  no  spurious  2Ax  or  2At  modes  due  to  the 
ability  of  the  GWCEFE  scheme  to  propagate  high  wave  number  energy.  Convergence 
properties  in  space  and  time  are  good  with  superconvergence  occurring  in  the  range 

Cr  =  0.5  -  1.5.  In  this  range,  more  accurate  solutions  are  obtained  using  larger  time 
steps. 

Application  of  ADCIRC-2DDI  to  the  English  Channel  and  Southern  North  Sea 

78.  The  accuracy  and  behavio’^al  characteristics  of  the  external  mode  solution 
have  been  tested  in  field  applications  including  (a)  tidal  and  hurricane  storm  surge 
simulations  in  the  Gulf  of  Mexico  (Westerink  et  al.,  in  review),  (b)  tidal  simulations 
in  the  English  Channel  and  Southern  North  Sea,  (c)  tidal  simulations  in  a  small 
coastal  inlet  (Luettich,  Birkhahn,  and  Westerink  1991)  and  (d)  tidal  simulations  in 
the  New  York  Bight.  The  English  Channel/Southern  North  Sea  system  is  probably 
the  best  documented  field  site  presently  in  existence  for  testing  a  long-wave, 
hydrodynamic  model.  Since  the  emphasis  of  this  report  is  on  the  development  and 
testing  of  the  various  components  of  ADCIRC,  the  results  of  applying  ADCIRC-2DDI 
to  the  English  Channel  and  Southern  North  Sea  are  presented  below. 

79.  In  the  mid-1980’s  considerable  effort  was  put  forth  to  establish  and  make 
readily  available  a  set  of  standard  grids,  boundary  conditions,  and  verification  data  for 
model  evaluation  in  the  English  Channel  and  Southern  North  Sea  (Werner  and  Lynch 
1988).  This  data  has  been  used  as  the  basis  for  modeling  studies  for  the  Tidal  Flow 
Forum  I  at  the  Conference  on  Finite  Elements  in  Water  Resources,  Lisbon,  Portugal, 
in  1986  and  for  the  Tidal  Flow  Forum  II  at  the  VII  International  Conference  on 
Computational  Methods  in  Water  Resources,  Cambridge,  MA  in  1988.  Two  collections 
of  scientific  papers  have  been  published  from  this  work  and  can  be  found  in  Advances 
in  Water  Resources.  Vol.  10,  No.  3  (1987)  and  Advances  in  Water  Resources.  Vol.  12, 
Nos.  3  and  4  (Dec  1989). 

80.  The  fully  nonlinear  version  of  ADCIRC-2DDI  was  applied  to  the  grid  and 
bathymetry  shown  in  Figure  10.  The  grid  consists  of  990  nodes  and  1,762  linear 
triangular  elements.  The  model  was  forced  by  specifying  11  harmonic  constituents  for 
elevation  (Oi,  K|,  M2,  N2,  S2,  K2,  MS4,  MN4,  M4,  Me,  2MSa)  along  the  two  open 
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modd  boundaries.  Wind  stress  and  tidal  potential  forcings  were  not  used  in  the 
modd  runs.  Modd  parameters  were  selected  to  match  those  used  by  previous 
investigators  to  allow  the  direct  comparison  of  model  results  with  field  data  and  with 
previously  published  modd  results.  The  following  parameter  values  were  used  in  the 
model:  u  =  0.0002s-S  At  =  360s,  Cf  =  0.002322,  f  =  0.0001 13341s ‘i,  and  Ehj  =  0.0. 
The  time  integration  coeffidents  in  the  GWCE  were  set  to  Oi  =  0.35,  02  =  0.30,  and 
03  =  0.35. 

81.  ADCIRC-2DDI  was  run  for  the  short-term  test  case  suggested  by  Werner 
and  Lynch  (1988)  covering  the  period  from  0  hr  on  15  March  1976  to  24  hr  on  17 
March  1976.  Werner  and  Lynch  (1987)  found  that  it  was  necessary  to  use  a 
minimum  bathymetric  depth  of  15  m  throughout  the  model  domain  to  avoid 
generating  negative  water  depths  during  their  simulations.  ADCIRC-2DDI  ran 
successfully  using  a  minimum  bathymetric  depth  as  smaU  as  10  m,  although  the 
simulated  results  were  highly  insensitive  to  this  change  at  the  19  locations  where 
observational  data  were  available  (see  Figures  11  and  12). 

82.  The  first  47  hr  10  min  of  the  simulation  were  used  as  a  transient  start-up 
period.  Figures  11  and  12  present  comparisons  between  modeled  time  series  and 
observed  time  series  of  free  surface  elevation  (at  11  stations)  and  depth-averaged 
current  speed  and  direction  (at  8  stations)  for  the  final  24  hr  50  min  of  the 
simulation.  (The  locations  of  the  elevation  and  velocity  stations  are  shown  in 
Figure  10.  The  observed  time  series  were  actually  reconstructed  from  11  primary  tidal 
constituents  at  each  station.  The  tidal  constituents  correspond  to  those  used  to  force 
the  model  open  boundaries  and  were  extracted  from  raw  time  series  at  each 
observation  station  using  harmonic  analysis.)  In  general,  the  model  does  a  good  job 
of  simulating  the  observed  results.  Some  of  the  differences  can  be  attributed  to  local 
topographic  and  bathymetric  effects  and  to  the  inherent  problems  associated  with 
representing  bottom  stress  in  a  depth-integrated  model.  Also,  Werner  and  Lynch 
(1989)  point  out  that  the  model  results  contain  harmonic  constituents,  generated  by 
nonlinear  interactions  within  the  domain,  that  are  not  included  in  the  reconstructed 
observed  time  series.  By  filtering  this  energy  out  of  the  model  results,  they  were  able 
to  reduce  the  average  difference  between  the  simulated  and  observed  surface  elevations 
by  approximately  40  percent.  The  worst  comparison  occurs  at  the  tidal  elevation 
station  at  Christchurch  and  is  at  least  partially  due  to  the  neglect  of  the  channel 
between  the  Isle  of  Wight  and  the  mainland  (located  approximately  25  km  east  of 
Christchurch)  in  the  model  grid. 

83.  ADCIRC-2DDI  was  also  run  for  the  long-term  test  case  suggested  by 
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observation  - - 15m  minimum  depth  . 10m  minimum  depth 


Figure  11.  (Sheet  2  of  3) 
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Figure  11.  (Sheet  3  of  3) 
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Figure  12.  Time  series  of  simulated  and  observed  tidal  speed  aud  direction  at  various  stations  (Continued) 


Werner  and  Lynch  (1988)  covering  the  190-day  period  starting  at  0  hr  on  15  March 
1976.  The  first  5  days  were  discarded  to  allow  for  start-up  transients  and  the 
remaining  185  days  were  harmonically  analyzed  using  the  least  squares  package  of 
Foreman  (1977).  The  amplitudes  and  phases  of  the  primary  surface  elevation 
constituents  from  the  ADCIRC  simulation,  from  a  simulation  by  Werner  and  Lynch 
(1989),  and  from  the  observed  time  series  at  the  11  elevation  stations  are  compared  in 
Table  3.  The  overall  comparison  between  model  results  and  observations  is  reasonable 
considering  no  effort  has  been  made  to  calibrate  the  model  by  adjusting  the  bottom 
friction  coefficient,  as  attempted  by  Baptista,  Westerink,  and  Turner  (1989).  Some  of 
the  largest  differences  in  phase  occur  at  stations  that  are  close  to  amphidromes.  This 
is  because  a  small  displacement  of  an  amphidrome’s  position  can  result  in  a  large 
change  in  the  nearby  phase  values.  Some  of  the  largest  relative  differences  in 
amplitude  (i.e.,  percent  difference  between  the  simulated  and  observed  amplitude) 
occur  in  the  Ma  constituents.  Bottom  friction  is  the  primary  nonlinear  generating 
mechanism  for  this  constituent,  suggesting  that  this  process  is  not  captured  very  well 
by  a  depth-integrated  model. 

84.  Figure  13  presents  co-tidal  charts  for  the  entire  domain  for  14  tidal 
constituents.  The  ADCIRC-2DDI  results  presented  in  Figure  13  and  Table  3  compare 
very  closely  with  those  of  Werner  and  Lynch  (1989).  This  is  expected  since  Werner 
and  Lynch  (1989)  used  a  depth-integrated,  finite  element,  GWCE-based  model  that  is 
similar  to  ADCIRC-2DDI.  The  minor  deviations  between  the  models  are  due  to 
ADCIRC’s  use  of  a  non-conservative  formulation  of  the  advective  terms  in  the  GWCE 
as  well  as  slight  differences  in  the  discretizations  of  several  of  the  terms.  The  close 
correspondence  between  the  model  results  provides  an  excellent  verification  of  the 
formulation  and  numerical  discretization  used  in  the  external  mode  of  ADCIRC. 
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a.  Oi  constituent 


Figure  13.  Co-tidai  charts  for  simnlated  constituents  (Sheet  1  of  7) 
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c.  Mj  constituent 


d.  Nj  constituent 


Figure  13.  (Sheet  2  of  7) 
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e.  Sj  constitueiit 


f.  Ka  constituent 


Figure  13.  (Sheet  3  of  7) 


g.  2MS2  constituent 


h.  2MN2  constituent 


Figure  13.  (Sheet  4  of  7) 


Ic.  MS4  constitneiit 


1.  Me  constituent 


Figure  13.  (Sheet  6  of  7) 


m.  2MSe  constituent 


n.  Residual  constituent 


Figure  13.  (Sheet  7  of  7) 
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PART  IV:  INTERNAL  MODE  SOLUTION 


Definition  and  ApplicabilitY  of  a  3DL  Model 

85.  As  discussed  in  Part  II,  mode  splitting  replaces  the  direct  solution  of  the 
three-dimensional  governing  equations  with  an  ”external  mode”  computation  for  free 
surface  elevation  (using  the  vertically  integrated  governing  equations)  and  an  "internal 
mode"  computation  for  the  vertical  profile  of  velocity.  It  was  noted  in  Part  II  that 
all  of  the  physics  contained  in  the  three-dimensional  governing  equations  are  included 
in  the  vertically  integrated  equations  if  the  bottom  stress  and  the  momentum 
dispersion  terms  are  specified  correctly.  The  simple  parameterizations  of  bottom  stress 
and  momentum  diSj,>ersion  in  terms  of  depth-averaged  velocity  (Equations  28  -  31)  are 
physically  correct  only  for  the  simplest  flows  (e.g.,  a  logarithmic  velocity  profile  over 
depth).  Mode  splitting  replaces  these  simple  parameterizations  with  the  internal  mode 
equations.  Therefore,  when  the  complete  internal  mode  equations  are  solved,  the 
bottom  stress  and  momentum  dispersion  used  in  the  vertically  integrated  equations  are 
(in  theory)  completely  consistent  with  the  three-dimensional  equations. 

86.  While  the  external  mode  equations  are  two-dimensional,  the  internal  mode 
equations  retain  the  spatial  variation  of  velocity  in  three  dimensions.  Considerable 
computational  savings  can  be  realized  if  the  advective  terms  and  the  horizontal 
momentum  diffusion  terms  are  dropped  in  the  internal  mode  computations  (Nihoul  and 
Djenidi  1987;  Davies  1988).  This  simplification  eliminates  all  horizontal  gradients  from 
the  internal  mode  equations,  thereby  reducing  them  to  one-dimensional  equations  in 
space  (over  the  vertical).  When  simplified  internal  mode  equations  are  solved,  the 
bottom  stress  and  momentum  dispersion  are  no  longer  completely  consistent  with  the 
three-dimensional  equations.  However,  these  approximations  should  be  physically 
correct  for  flows  in  which  the  vertical  distribution  of  momentum  at  each  horizontal 
grid  point  is  determined  by  a  local  balance  between  the  surface  and  bottom  stresses, 
vertical  momentum  diffusion,  the  Coriolis  force,  and  the  local  inertia.  (Clearly,  this 
should  encompass  a  much  wider  range  of  flows  than  parameterizations  solely  in  terms 
of  the  depth-averaged  velocity.)  The  required  balance  will  exist  when  the  rate  of 
vertical  momentum  transport  is  much  greater  than  the  rate  of  horizontal  momentum 
transport.  Assuming  horizontal  momentum  transport  is  dominated  by  advection,  the 
rate  of  vertical  momentum  transport  will  be  much  greater  than  the  rate  of  horizontal 
momentum  transport  in  the  three-dimensional  governing  equations  if 
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»  u 


da 

m 


Scaling  this  yields 


EvcUc 


» 


Ug 


Lc  ^  ^  hcUc 
or  TT"  ~T? — 

He  VC 


where  Eye,  Uc,  he,  and  Lc  are  a  characteristic  vertical  eddy  viscosity,  horizontal 
velocity,  water  depth,  and  horizontal  length  scale,  respectively.  Dimensional  arguments 
suggest  Eve  »  ^cUc  where  ^  is  a  constant  whose  value  for  tidal  and  wind-driven 
flows  typically  ranges  from  10‘3  to  10*2  (Bowden,  Fairmairn,  and  Huges  1959;  Csanady 
1976;  Fischer  et  al.  1979;  Davies  1985).  Therefore,  the  simplified  internal  mode 
equations  should  be  an  accurate  approximation  to  the  full  internal  mode  equations 
provided 

^  »  100  -  1,000 

Since  coastal  and  shelf  waters  are  usually  characterized  by  large  length-to-depth 
scales,  a  model  based  on  the  simplified  internal  mode  equations  should  be  widely 
applicable  in  these  waters. 

87.  The  model  based  on  the  simplified  internal  mode  equations  will  be  called  a 
three-dimensional  local  (SDL)  model.  This  name  emphasizes  the  fact  that  the 
simplified  internal  mode  equations  give  values  of  bottom  stress  and  momentum 
dispersion  for  the  two-dimensional  (external  mode)  equations  that  are  not  fully 
consistent  with  three-dimensional  equations,  but  rather  are  based  on  a  local 
approximation  of  the  three-dimensional  equations. 


Rationale  for  the  DSS  Technique 


88-  Despite  the  savings  gained  by  simplifying  the  internal  mode  equations  in 
the  SDL  model,  the  internal  mode  equations  are  difficult  to  solve  numerically  because 
of  the  high  velocity  gradients  that  characterize  the  water  column  near  the  bottom  and 
surface  boundaries  and  across  strong  density  changes.  Existing  state-of-the-art 
circulation  models  use  velocity  as  a  dependent  variable  and  therefore  require  a  fine 
numerical  discretization  to  resolve  regions  of  rapid  velocity  change.  Davies  (1991)  and 
Davies  and  Jones  (1991)  have  examined  the  computational  effort  required  to  resolve  a 
bottom  boundary  layer  using  a  one-dimensional  model  through  the  vertical  solved  with 
finite  differences  and  several  coordinate  transformation/grid  stretching  schemes.  For 
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tidal  flows  having  an  eddy  viscosity  that  is  constant  over  the  upper  80  percent  of  the 
water  column  and  that  varies  linearly  with  distance  from  the  bed  over  the  bottom  20 
percent,  Davies  (1991)  found  that  it  was  necessary  to  use  a  logarithmic  or  log-linear 
coordinate  transformation  and  at  least  20  grid  cells  to  obtain  convergence  of  the 
velocity  sclation.  When  the  eddy  viscosity  was  determined  using  a  level  2-1/2 
turbulent  closure,  the  most  efficient  solution  was  found  to  require  a  log^inear 
coordinate  transformation  and  50  -  100  grid  cells  over  the  vertical  for  both  a 
turbulent  kinetic  energy  transport  equation  and  the  momentum  equation. 

89.  Practical  geophysical  flows  often  have  two  or  more  regions  containing  sharp 
velocity  gradients  over  the  vertical.  Because  of  the  computational  overhead  in  time 
and  memory  required  to  resolve  these  features,  existing  multi-<limensionai  circulation 
models  almost  always  omit  the  near  bottom  region  and  use  a  slip  boundary  condition 
that  expresses  bottom  stress  as  a  quadratic  function  of  near  bottom  velocity.  This 
assumption  is  physically  correct  only  when  the  velocity  profile  below  the  lowest  grid 
point  is  logarithmic.  An  accurate  treatment  of  surface  and/or  internal  boundary 
layers  requires  a  fine  grid  in  the  regions  of  these  layers.  In  many  cases  the  required 
computational  overhead  makes  it  impractical  to  resolve  these  features  in  multi¬ 
dimensional  computations.  A  survey  of  the  recent  literature  suggests  that  only  rarely 
have  more  than  ~  20  grid  cells  been  used  over  the  vertical  in  three-dimensional 
engineering  or  geophysical  model  applications.  For  example,  Oey,  Mellor,  and  Hires 
(1985^  used  11  grid  cells  over  the  vertical  in  their  model  of  the  Hudson-Raritan 
Estuary.  Clearly,  such  models  have  limited  ability  to  resolve  even  one  significantly 
sheared  velocity  gradient  region.  (Note:  Davies  and  Jones  (1990)  have  recently 
published  results  from  a  three-dimensional  model  of  the  northern  European  continental 
shelf  using  45  grid  cells  over  the  depth.  However,  this  model  uses  a  coarse  horizontal 
grid  and  omits  the  advective  terms  in  both  the  internal  and  the  external  mode- 
governing  equations.) 

90.  It  is  well-established  from  laboratory  and  field  experiments,  theoretical 
arguments,  and  conventional  one-dimensional  models  that  the  time-averaged  vertical 
shear  stress  varies  rather  smoothly  through  the  water  column,  particularly  near 
boundaries.  Therefore,  it  should  be  possible  to  use  a  relatively  coarse  vertical 
discretization  to  solve  numerically  for  the  vertical  shear  stress,  even  in  boundary 
layers.  A  novel  technique  has  been  developed  that  allows  the  vertical  shear  stress  to 
be  used  in  place  of  velocity  as  the  dependent  variable  in  the  internal  mode  equations. 
Applications  of  the  DSS  technique  using  linearized  equations  of  motion  (discussed  in 
detail  below)  have  shown  that  it  provides  a  highly  efficient  means  of  solving  the 
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internal  mode  equations.  This  technique  promises  to  be  invaluable  for  modeling 
coastal  and  shdf  circulation  in  which  the  bottom  and  surface  boundary  layers  comprise 
a  significant  portion  of  the  water  column  and  for  modeling  processes  that  are  critically 
dependent  on  boundary  layer  physics  such  as  wave-current  interaction,  sediment 
transport,  oil  spill  movement,  ice  floe  movement,  energy  dissipation,  physical-biological 
couplings,  etc. 


Devdopment  and  Testing  of  DSS  Method  No.  1 

91.  Internal  mode  equations  can  be  generated  by  subtracting  the  vertically 
integrated  equations  firom  the  three-dimensional  equations  (Wang  1982;  Sheng  1983; 
Davies  1985).  Using  the  three-dimensional  equations  in  the  a  coordinate  system 
(Equations  19  -  21),  the  non-conservative  vertically  integrated  momentum  equations 
(Equations  25  and  26),  assuming  a  constant  density  fluid,  and  neglecting  advection 
and  horizontal  momentum  diffusion  terms,  the  resulting  internal  mode  equations  are 

92.  Using  the  eddy  viscosity  rdationships  (Equation  34)  to  express  Tzx  and  r^y 
in  terms  of  velocity  and  either  the  slip  or  the  no-slip  boundary  condition  (Equation 
10)  at  the  bottom.  Equations  114  and  115  can  be  cast  entirely  in  terms  of  velocity. 
Numerical  solutions  can  then  be  sought  for  the  dependent  variables,  u  and  v.  This  is 
the  standard  velocity  solution  (VS)  approach. 

93.  Alternatively,  Equation  34  can  be  inverted  to  obtain  expressions  for 
velocity  in  terms  of  stress 

»  =  -  u  + 

b 

V  =  Vb  -  V  + 

b 

In  Equation  116  the  definitions  of  u  and  v  have  been  used  and  nonzero  slip  velocitira 
Ub  and  Vb  have  been  included  for  generality.  Relating  Ub  and  vb  to  the  bottom 
stresses,  rb*  and  Tby,  via  the  slip  conditions 
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rbx/Po  =  k  Ub  =  k(ab  +  U) 
Tby/po  =  k  Vb  =  k(vb  +  V) 


(117a) 

(117b) 


Equation  116  can  be  written  as 

=  Ibx  _  u  +  I  da 
Po^  (a— bj  J  PqEv 

b 

'  =  ^  I  ^ 


(118a) 

(118b) 


(For  a  no-slip  boundary  condition,  the  terms  and  do  not  appear  in  Equation 

118.  The  no-slip  condition  is  approached  as  k-^oo.)  Substituting  Equation  118  into 
Equations  114  and  115  gives: 


]  = 


jry  Ts 


FT 


(119) 


d  H 


5T 1  J  ^ 

b  b 

“  +  ^y]  =  It  ^  ^  ^^20) 


Equations  119  and  120  have  tzx  and  Tzy  as  dependent  variables  and  will  be  called  the 
DSS‘  internal  mode  equations.  (The  superscript  1  is  used  to  identify  DSS  method 
No.  1.)  These  equations  are  forced  by  the  external  mode  solution  (U,  V,  and  H)  and 
the  applied  surface  stress. 

94.  Equations  119  and  120  contain  both  integral  and  differential  terms; 
therefore,  they  are  well-suited  for  a  spatial  discretization  in  which  Tzx  and  Tzy  are 
expressed  in  terms  of  assumed  shape  functions  such  as  the  spectral  or  finite  element 
methods.  Depending  on  the  choice  of  the  shape  functions  and  the  functional  variation 
of  Ey  over  the  depth,  the  velocity  profile  can  be  recovered  from  the  stress  profile  by 
solving  Equation  118  in  closed  form.  This  is  an  important  convenience  because  it 
avoids  the  troublesome  operation  of  numerically  integrating  the  near^ogarithmic 
singularity  that  occurs  in  Equation  118  when  Ey  varies  with  distance  from  a 
boundary.  The  restrictions  that  a  closed-form  solution  for  Equation  118  impose  on 
T’zx)  Tzy,  and  Ey  are  not  severe.  For  example,  Tzx,  Tzy,  and  Ey  may  be  expressed  in 


terms  of  polynomials  that  span  the  vertical  globally  or  in  a  piecewise  manner. 
Polynomial  variations  of  Tzx  and  r^y  are  consistent  with  either  the  spectral  method  or 
the  finite  element  method;  for  most  practical  problems,  Ey  can  be  approximated  as 
piecewise  linear  over  the  vertical  (Fumes  1983;  Davies  1987;  Chu,  Liou,  and  Flenniken 
1989;  Jenter  and  Madsen  1989). 

95.  The  effectiveness  of  the  DSS>  technique  is  evaluated  using  a  simple  test 
case  consisting  of  flow  generated  by  a  specified  surface  stress  aligned  in  the  x-direction 
in  a  wide,  straight  channel  of  constant  depth  with  no  Coriolis  force.  An  analytical 
solution  can  be  found  for  the  linear  version  of  this  problem  and  provides  a  benchmark 
for  the  numerical  solutions.  For  convenience,  the  linear  governing  equations  are 
repeated  below: 


dU 


(121) 

(122) 

(123) 

(124) 


Equations  121  and  122  are  the  depth-integrated  (external  mode)  continuity  and 
momentum  equations;  Equation  123  is  the  VS  internal  mode  equation;  Equation  124  is 
the  DSS‘  internal  mode  equation.  Since  there  is  no  motion  in  the  y-direction,  the 
y-direction  equations  and  the  subscript  "x"  in  the  stress  terms  have  been  dropped. 

96.  The  Galerkin-spectral  method,  with  shape  functions  consisting  of  Legendre 
polynomials  (LPs)  over  the  interval  -1  <  cr  <  1  is  used  to  discretize  the  VS  and  the 
DSS*  internal  mode  equations.  The  order  LP  is  denoted  Lb  and  can  be  computed 
from  the  recursion  formulas 
Lo(£r)  =  1 
Li{a)  =  ff 

=  [^]  fir  -  Lr-i 

The  first  eight  LPs  are  shown  in  Figure  14.  Other  properties  of  LPs  of  note  are 
Lrs  =  Lr(l)  =  1 
Lrb  =  Lr('l)  =  (-I)-- 

d(r  =  2 
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Polynomials 


I^Lr(o^)  d<r  =  0  for  ail  r  >  1 


It  has  been  shown  for  wind-driven  circulation  that  velocity  solutions  using  Legendre 
and  Chebyshev  polynomials  yield  results  of  virtually  identical  accuracy,  that  these  are 
highly  superior  to  velocity  solutions  obtained  using  expansions  of  trigonometric 
functions,  and  that  these  are  more  accurate  than  velocity  solutions  computed  with  a 
second-order  finite  difference  scheme  having  the  same  number  of  degrees  of  freedom 
(Davies  and  Owen  1979;  Davies  and  Stephens  1983).  For  further  information  on  the 
use  of  spectral  methods  in  three-dimensional  circulation  models,  the  interested  reader 
is  referred  to  an  excellent  review  by  Davies  (1987). 

97.  The  Galerkin-spectral  discretization  for  the  VS  internal  mode  equation  is 
obtained  by  multipipng  Equation  123  by  the  weighting  fuiiCtion  and  integrating 
from  -1  to  1,  i.e., 

4|‘t.  u  d.  -  I  (l,  y.^)  d,  =  -  J  [a  -  a]  |‘i.  d,  (125) 

-1  -1  -1 

Integrating  the  second  term  in  Equation  125  by  parts 


Ih 

Po 


1 


Tz 

Po  ^ 


(126) 


and  substituting  this  into  Equation  125  yields; 

If  j  I,  u  d(r  +  I  j  g  d<7  =  -  J  j  L.  der 

-i  -1  -1 

[ 


-I 


T  Ts. 


Using  the  definition  of  the  LP,  Equation  127  simplifies  to 


0  =  0 
1 


If  I  L„  u  d(r  +  I  I  ^  |§a  da  -  ^  -  La,b 

-1  -1 


Ik 

Po 


(127) 

m  =  0 

(128) 

m  >  1 

(129) 

Since  Lo{a)  =  1,  the  operation  that  generates  Equation  128  is  equivalent  to 
integrating  Equation  123  over  the  depth  when  m  =  0.  The  identity  in  Equation  128 
occurs  because  Equation  123,  by  definition,  has  no  depth-averaged  component. 
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98.  The  final  steps  in  applying  the  Galerkin-spectial  discretization  to  the  VS 
are  to  substitute  Equation  34  for  Tz/po  in  Equation  129  (noting  that  ^  and  to 

expand  u  as  a  series  of  LPs  with  time-varying  coefficients,  ^n(t),  i.e., 

R 

u(ff,t)  =  J^(t)  Ln((r)  (130) 

n*t 

I 1  I  p^]  ®  ^ 

n»l  -1  n»l  -I 


Because  |'Z/n  d(r  =  0  for  n  >  1,  the  necessary  condition  |'u  da  =  0  is  identically 

satisfied  by  the  spectral  solution  by  using  only  the  n  >  1  LP.  The  solution  of 
Equation  131  requires  a  bottom  boundary  condition  (Equation  117).  After  expanding 
u,  this  becomes 


R 

I  ^nb  =  -  U  + 

n  =  t 


(132) 


99.  The  Galerkin-spectral  discretization  for  the  DSS‘  is  obtained  by 
multiplying  Equation  124  by  the  weighting  function  Ln{a)  and  integrating  from  -1  to 
1,  i.e., 


lit  [f  I  ^  [t  It^^)  ■  fe]  I  da  -  I  j  L.  |5(^)  da  - 
-1  -1  -1  -1 


.It  -  Efc]  1 ' 

-1 


Lm  da 


Integrating  the  stress  derivative  by  parts  changes  Equation  133  to 


[f  ^  1  ‘‘‘'‘‘"J  +  [e  Ie<S) 


h  d 

J  ^IJ  J 

-t  -1 


Using  the  definition  of  the  LP,  Equation  134  simplifies  to 


(133) 


-  E  ^  ^  - 1  ^  -  [It  -  ffe]  I  ^ 

-1  -1 
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0  =  0 


m  =  0 


(135) 


jltUi- +  +  =  036) 


-I  -i 


The  id^tity  in  Equation  135  occurs  because  Equation  136  has  no  depth-averaged 
component. 

100.  The  final  step  in  applying  the  Galerkin-spectral  method  for  the  DSS*  is 
to  expand  rz/po  as  a  series  of  LPs  with  time-varying  coefficients,  an(t),  i.e., 


^  ^  I «»(‘)  i«(^) 


»  1 


»  1 


I"  I  l?^|  ^  ^  I  J  ^nbfims  ]  - 


(137) 


m  >  1  (138) 


n  =  0  '-i 


n  =  0  -I 


101.  The  bottom  boundary  condition  was  introduced  into  Equation  118  and 
subsequently  into  Equation  138.  Therefore  it  does  not  generate  an  extra  equation,  as 
was  the  case  for  the  VS.  However,  the  stress  expansion,  Equation  137,  does  not 
automatically  satisfy  the  condition  that  ^  u  da  =  0.  Rather  this  must  be  enforced 

explicitly.  Using  Equation  118  and  the  definition  of  u,  this  requirement  generates  the 
additional  equation 


h 


da  =  0 


(139) 


Substituting  the  expansion  for  rz/po  into  Equation  139  yields 


(140) 


102.  The  relative  merit  of  the  DSS‘  versus  the  VS  was  evaluated  by 
comparing  solutions  computed  numerically  with  analytical  solutions  for  the  problem  of 
wind-driven  circulation  in  a  closed,  rectangular  channel  aligned  with  the  x-axis  and 
having  a  constant  bathymetric  depth.  This  was  done  for  a  steady-state  case,  for  a 
periodically  varying  wind  stress,  and  for  an  instantaneously  imposed  wind  stress. 

103.  In  each  test  case,  Ey  was  assumed  to  be  linear  over  the  depth  as 
expressed  by 

Ev(a)  =  Ezo(a+l4-ao)  (141) 
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where  Cq  =  2zo/h  is  the  dimensionless  roughness  height.  It  is  well  known  &om 
theoretical,  laboratory,  and  field  experiments  that  an  eddy  viscosity  that  increases 
linearly  with  distance  from  a  solid  boundary  realistically  reproduces  the  physics  of  the 
boundary  layer  near  the  boundary  (Monin  and  Yaglom  1971;  Schlichting  1979;  Grant 
and  Madsen  1986).  Despite  the  fact  that  this  does  not  hold  over  the  entire  depth, 
(e.g.,  it  has  been  suggested  that  Ey  should  also  increase  linearly  with  distance  below 
the  free  surface  (Jenter  and  Madsen  1989)),  Equation  141  is  used  here  because  it 
generates  a  realistic  bottom  boundary  layer  and  because  it  simplifies  the  analyses  of 
model  results  by  introducing  only  two  parameters,  Ego  and  Co,  into  the  problem.  As 
is  shown  below,  the  presence  of  a  velocity  gradient  region  at  the  bottom  is  sufficient 
to  illustrate  the  advantage  of  the  DSS  over  the  VS.  In  fact,  the  use  of  an  eddy 
viscosity  that  does  not  also  give  a  boundary  layer  at  the  free  surface  is  a  considerable 
advantage  for  the  VS,  since  it  eliminates  the  additional  need  to  reproduce  velocity 
gradients  there. 

104.  Assuming  reasonable  ranges  for  Zq  of  0.1  to  10  cm,  and  for  h  of  1  to 
100  m,  suggests  values  of  ao  ~  10**  to  lO'*.  (The  combination  of  Zq  =  10  cm  and 
h  =  1  m,  which  gives  Cq  ~  10**,  is  not  considered  realistic  since  Zq  is  typically 

3  to  10  percent  of  the  physical  roughness  height.  In  this  case  the  physical  roughness 

would  occupy  the  entire  depth.)  Assuming  the  slope  of  the  variation  of  Ey  with  z 
*  *  * 

scales  with  Ub,  (Ub  =  V  n)/po  ),  then  Ego  ~  Ubh.  If  Ub  varies  over  the  range  0.1  to 
10  cm/s,  Ego  ~  10*5  to  10*  m^/s. 

105.  Equations  131,  132,  138,  and  140  show  that  the  VS  and  DSS*  require  the 
specification  of  Ta/po  (which  is  the  input  forcing)  and  U.  To  eliminate  the  possibility 
that  errors  in  the  solution  for  U  might  affect  the  comparisons,  U  was  obtained  for 
each  test  case  from  an  analytical  solution  of  Equations  121  -  123.  As  a  result,  errors 
in  the  VS  and  DSS*  over  the  vertical  do  not  feed  back  into  the  solution  for  U  as 
they  would  if  the  complete  problem  was  solved  numerically. 

106.  In  all  of  the  results  presented  below,  bottom  stresses  are  obtained  from 
the  VS  by  using  computed  bottom  slip  velocities  and  the  linear  slip  boundary 
condition  (Equation  117).  Comparisons  indicated  that  this  method  gave  more  accurate 
values  of  bottom  stress  than  those  obtained  by  evaluating  Equation  34  at  a  =  -1. 

(A  similar  conclusion  was  reached  by  Gresho,  Lee,  and  Sani  (1987).)  Velocities  are 
obtained  from  the  DSS*  by  solving  Equation  118  analytically  using  the  computed  stress 
profiles. 
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(142) 


108.  The  VS  and  DSS^  are  obtained  from  Equations  131,  132,  138,  and  140  by 
dropping  the  time  derivatives,  setting  U  =  0,  and  considering  all  other  terms  to  be 
constant  in  time. 

109.  Figure  15  presents  a  comparison  of  vertical  profiles  of  horizontal  velocity 
for  several  combinations  of  K  and  Co  computed  from  the  analytical  solution,  the  DSS* 
using  2  LPs  and  the  VS  using  various  numbers  of  LPs.  Equation  143  indicates  that 

Ok 

the  analytical  solution  for  stress  varies  linearly  over  the  depth,  regardless  of  the  form 
of  Ey.  This  solution  can  be  represented  exactly  by  the  DSS‘  using  only  the  n  =  0 
and  n  =  i  LP;  therefore  the  DSS'  and  the  analytical  solution  in  Figure  15  are 
identical.  Equation  146  indicates  that  the  analytical  solution  for  velocity  has  a 
logarithmic  variation  over  the  depth  and  consequently  a  potentially  sharp  gradient 
region  near  the  bottom.  In  Figure  15a  the  combination  of  a  small  K  (large  amount 
of  slip)  and  a  large  <7©  minimizes  the  gradient  region.  Over  most  of  the  depth  the 
velocity  profile  is  nearly  linear  and  therefore  closely  reproduced  using  a  VS  with  2  LP. 
However,  approximately  5  LPs  are  required  to  capture  the  mild  velocity  gradient  near 
the  bed.  In  Figure  15b,  the  same  K  is  used  with  <7o  reduced  by  two  orders  of 
magnitude.  This  has  the  effect  of  pushing  the  gradient  region  closer  to  the  bottom 
(i.e.,  it  is  equivalent  to  increasing  the  depth  by  a  factor  of  100  for  the  same 
roughness)  and  therefore  steepening  the  velocity  gradient.  Because  the  velocity  profile 


100 


Solution  -  DSS  2  LP  (o)  1  I  — Exoct  Scrfution  -  DSS  2 


Vertical  profiles  of  horizontal  velocity  for  the  steady  state  test  case  for  the  VS  and  DSS^ 


is  nearly  linear  over  much  of  the  depth,  it  is  reproduced  well  by  the  VS  with  2  LPs. 
However,  near  the  bed,  approximately  10  LPs  are  required  for  the  V3  to  capture  the 
gradient  region.  As  discussed  below,  this  results  in  a  poor  prediction  of  bottom  stress. 

110.  In  Figures  15c  and  15d,  a  high  value  of  K  is  used,  resulting  in  essentially 
no  slip  at  the  bottom.  For  large  <To  (Figure  15c)  a  velocity  expansion  of  10  or  more 
LPs  is  required  to  reproduce  this  profile.  Reducing  ao  by  two  orders  of  magnitude 
(Figure  15d)  sharpens  the  profile  further,  and  approximately  20  LPs  are  required  to 
capture  the  velocity  profile  away  from  the  boundary.  Many  more  are  required  to 
represent  the  gradient  region  near  the  boundary. 

111.  As  noted  above,  an  important  reason  for  using  a  three-dimensional  model 
in  place  of  a  two-dimensional  model  is  the  former’s  improved  representation  of  the 
bottom  stress.  However,  since  stress  is  proportional  to  the  velocity  gradient 
(Equation  34),  or  the  bottom  velocity  (Equation  117),  the  bottom  stress  may  still  be 
represented  poorly  if  the  gradient  region  near  the  bottom  is  not  resolved  properly.  To 
illustrate  this  problem,  a  comparison  was  made  between  the  analytical  bottom  stress 
and  computed  bottom  stresses  from  the  DSS‘  and  the  VS  over  the  practical  range  of 
K  and  ao-  The  DSS*  reproduces  bottom  stress  exactly  using  2  LPs.  On  the  other 
hand,  Table  4  presents  a  summary  of  the  number  of  LPs  requited  for  the  computed 
bottom  stress  using  the  VS  to  come  within  10  percent  of  the  analytical  bottom  stress 
as  a  function  of  K  and  ao-  Clearly,  it  is  computationally  practical  to  use  the  VS 
only  for  large  roughnesses  and  large  amounts  of  slip,  both  of  which  tend  to  minimize 
the  velocity  gradient  at  the  bottom. 

112.  Although  quite  simple,  the  steady-state  case  demonstrates  the  relative 
ease  with  which  a  DSS  can  resolve  a  realistic  boundary  layer  (i.e.,  no  bottom  slip  and 
a  linearly  varying  eddy  viscosity)  in  a  hydrodynamic  model  that  explicitly  includes  the 
vertical  dimension.  In  the  following  examples  we  evaluate  how  this  highly  desirable 
capability  is  affected  by  unsteady  conditions.  Only  the  no-slip  case  (K  =  1,000)  is 
considered. 

113.  If  a  periodic  surface  stress  is  assumed  of  the  form  rs(t)//7o= 

(where  u  is  the  forcing  frequency  and  i  5  solutions  can  be  sought  to 

Equations  121  -  123  that  have  the  form  U(t)  =  Ueiu>*i,  u(£r,t)  =  u(a)ei‘^^  rb(t)/po= 
{'Th/ stnd  r/(t)  =  (Note:  Tg/po,  U,  u(«r),  rb/po,  and  are  all  complex 

variables;  therefore  they  may  be  out  of  phase  with  each  other.)  Substituting  these 
into  Equations  121  -  123  transforms  the  linear  hydrodynamic  equations  into 
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Table  4 


Stftadv-State  Bottom  Stresses  Computed  Using  Velocity  Expansiong 


7"b(  anal -llx  c 

K 

#LP  J 

■Hjc  anal ) 

10-2 

10-1 

3 

0.100 

10-2 

10® 

8 

0.091 

10-2 

101 

g 

0.099 

10-2 

102 

10 

0.078 

10-2 

103 

10 

0.078 

10-3 

10-1 

8 

0.096 

10-3 

100 

21 

0.098 

10-3 

101 

24 

0.095 

10-3 

102 

24 

0.098 

10-3 

103 

24 

0.099 

10-4 

10-1 

22 

0.100 

10-4 

lO® 

<40 

0.192* 

10-4 

101 

<40 

0.242* 

10-4 

102 

<40 

0.249* 

10-4 

103 

<40 

0.249* 

10-5 

10-1 

<40 

0.174* 

10-5 

10® 

<40 

0.476* 

10-5 

101 

<40 

0.602* 

10-5 

102 

<40 

0.619* 

10-5 

103 

0.620* 

*  This  is  the  minimum  difference  obtained  using  no  more  than  40  Legend*o  polynomials. 
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(147) 

(148) 


iu}fj  +  h  =  0 


1 

•wU  =  -g  (r,  -  Tb) 

iu4  -  |-  f(»-yv 


35 


(149) 


114.  The  procedure  used  to  solve  Equations  147  -  149  analytically,  together 
with  the  linear  slip  boundary  condition,  has  been  presented  previously  (Lynch  and 
Officer  1985;  Lynch  and  Werner  1987)  and  is  not  repeated  here.  Rather,  the  solutions 
are  given  without  derivation  in  Table  5  (Equations  150  -  163). 

115.  Spectral  approximations  for  the  periodic  case  are  generated  by  expressing 

/9n(t)  =  and  an(t)  =  and  substituting  these  as  well  as  the  periodic  forms 

of  u(t),  rz(<7,t)/po,  Ts(t)/po,  nj(t)/po  and  U(t)  into  Equations  130  -  132,  137,  138,  and 
140.  The  resulting  equations  for  the  VS  ate 


X 


u(cr) 

n  9  f 

ri  Ln( 

(164) 

X 

n  *  1 

.1 

X 

1 

n»l 

J 

LjuL-, 

■i 

n  d<r 

+ 

n»l 

J  W 

-l 

dLju  J 

W 

1 

«|  o 

II 

?  i^.b 

Po  j 

m  >  1  (165) 

H 

i0, 

nsl 

f'nb 

=  - 

U  + 

Th 

^0 

(166) 

and  for  the  DSS‘  are 

Po 

X 

=  ^  ®n  Ljj 

n  9  A 

(167) 

X 

M  —  V 

1 

H  1 

h^  V  . 

^  2 
n»0 

i- 

-1 

f  in 

1  e; 
- 1 

dffdff  + 

1 

I  ®n[f 

rt*0  -1 

^  d<T 

d 

+  LnbLmbj 

_  Is 

Po 

m  >  1  (168) 

X 

I 

IB  A 

■‘rt 

J  ^ 

-  dtrderj 

U 

-  ¥ 

(169) 

n*0  -1-1 


116.  The  periodic  solution  depends  on  the  dimensionless  parameters  K  and  Cq 
(as  found  for  the  steady-state  solution),  a  dimensionless  channel  length  L' ,  a 
dimensionless  frequency  fi,  and  the  dimensionless  position  in  the  channel  x/L.  L'  is 
the  ratio  of  the  channel  length,  L,  to  the  wave  length  of  a  shallow-water  wave  having 
period  u)  (Equation  154).  fl  is  the  ratio  of  the  time  scale  lor  momentum  to  be 
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Table  5 

Analytical  Solution  for  the  Periodic  Test  Case* 


UqEzo 

■Ai(a)  -  Ai' 

,  1 

"A 2(a)  -  ki 

hTg/po 

B 

2(2+(To) 

L  ^ 

UEzO  _  UqEzO 
h  TsIpq  hrg/pol 


1  + 


2(2+'^)  F- 


(l-exp(-AL'))exp(AL'^)  -  (1-exp (AL'))exp(-AL'^) 


UqEzo  _  i(lj-72) 
bTg/po  n 

L'  =  — 


exp(AL')  -  exp{-AL') 


vP" 


tiizo 


A  =  V  i7i  -  1 
n(Ai+B) 


72  = 


Ai(-1  )A2  _  A2(-l)'j 


B(Ai+B) 


B 


ao_ 

2+0-0 


Ai(o)  =  /xi(cr)/i2(l)  -  /X2(«^)Ai(1) 

A2(o^)  =  /ii(a)^^2(-l)  -  A2(-l)j  -  M2(o^)|^/ii(-l)  -  /fiC-l)! 

S  -  Ml(l)|^y^2(-1)  -  /*2(-l)j  -  /t2(l)^/^l(-l)  - 

fii{a)  =  ber^[n(a+l+Oo)]‘^2j  +  i  bei^[n(o-+l+Oo)]<^*j 
^l2{a)  =  ker  ^[n(o+l+cro)]^2j  +  i  kei^[n(o-+l+o-o)]’/*j 


(150) 

(151) 

(152) 

(153) 

(154) 

(155) 

(156) 

(157) 

(158) 

(159) 

(160) 

(161) 

(162) 

(163) 


*ber,  bei,  ker,  kd  are  zeroth  order  Kelvin  functions,  an  overdot  (  )  =  djda, 
an  overbar  (  )  =  j  /j  da 
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transported  through  the  water  depth,  h^/Ezo,  and  the  forcing  time  scale,  l/w, 

Equation  155.  Assuming  ranges  for  u  of  lO'S  sec-i  to  lO'®  sec'i,  L  of  1  km  to 
10®  km,  and  h  and  Ub  as  given  previously,  suggests  L'  ~  10'®  to  10®  and  fl  ~  10 to 
10®.  In  all  cases  results  are  presented  for  x/L  =  0.5,  as  these  are  representative  of 
the  behavior  throughout  the  rest  of  the  channel. 

117.  Figures  16  and  17  present  magnitude  and  phase  portraits  of  the  velocity 
structure  for  K  =  1,000,  L'  =  1,  and  four  combinations  of  ao  and  n.  For  the  case 
fl  =  10*1,  momentum  is  transported  through  the  depth  in  only  a  fraction  of  the 
forcing  period.  Figures  16a  and  16b  and  17a  and  17b  show  that  the  velocity 
magnitude  and  phase  obtained  from  the  DSS^  using  2  LPs  are  virtually  identical  to 
the  analytical  solution;  therefore,  the  stress  variation  is  very  close  to  linear  over  the 
depth.  This  linear  stress  variation  suggests  that  the  momentum  balance  over  the 
depth  is  nearly  at  steady  state  and  is  consistent  with  the  low  value  of  fl.  Since 
steady  state  is  approached  as  fl  -*  0,  the  DSS‘  using  two  LPs  gives  a  highly  accurate 
solution  for  fl  <  10*i  as  well.  The  VS  is  able  to  capture  the  phase  change  through 
the  water  column  with  a  comparable  number  of  LP  to  the  DSS*.  However,  as  was 
the  case  at  steady  state,  for  ao  =  10'^,  approximately  10  LPs  are  required  to 
reproduce  the  velocity  magnitude  with  an  accuracy  comparable  to  the  DSS*  using 

2  LPs.  For  (To  =  10'^,  more  than  20  LPs  are  required. 

118.  For  the  case  fl  =  10,  the  vertical  momentum  balance  is  no  longer  near 
steady  state;  consequently  the  DSS‘  requires  more  than  2  LPs  to  capture  the  vertical 
stress  variation.  Figures  16c  and  16d  and  17c  and  17d  suggest  that  approximately 

4  LPs  may  be  needed  by  the  DSS'.  The  VS,  however,  requires  at  least  10  LPs  for 
ao  =  10**,  and  more  than  20  LPs  for  ao  =  10*^. 

119.  Figures  18  and  19  compare  the  amplitude  and  phase  behavior  of  the 
analytical  solution  for  bottom  stress  with  solutions  obtained  using  the  DSS*  and  VS. 
These  runs  were  made  using  a  single  value  of  ao  =  10'*,  but  varying  fl  and  L'.  The 
10^  change  in  L'  has  minimal  effect  in  these  pictures,  indicating  that  the  number  of 
LPs  required  for  the  DSS*  or  the  VS  to  converge  to  the  analytical  solution  is  only 
very  weakly  dependent  on  L'.  For  fl  <  1,  the  DSS*  with  2  LPs  is  nearly  identical  to 
the  analytical  solution,  while  for  larger  fl  the  number  of  LPs  required  by  the  DSS* 
increases  to  as  many  as  7  for  fl  =  10®.  Considering  the  fact  that  comparable  results 
using  the  VS  require  the  use  of  more  than  20  LPs,  the  DSS’  is  computationally  quite 
superior  to  the  VS  for  all  fl. 

120.  Although  the  Coriolis  force  was  omitted  from  these  test  cases,  the  results 
can  be  used  to  infer  whether  a  DSS  will  be  equally  effective  when  the  Coriolis  force  is 
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Pirure  16.  Vertical  profiles  of  the  horizontal  velocity  magnitude  for  the  periodic  test  case  for  the 
^  VS  and  DSSi 
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Vertical  profiles  of  the  horizontal  velocity  phase  for  the  periodic  test  case  for  the  VS  and  DSS‘ 
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Figure  19.  Ratio  of  the  spectral  to 
VS  and  DSS^ 


included.  The  counterparts  to  Equation  123  for  the  case  in  which  Coriolis  is  included 
are 


iwii  -  fv  =  _  r*,)  (170) 

fii  =  ((a-b^Ev  ^  ^  ^171) 

It  has  been  shown  (Lynch  and  Officer  1985)  that  the  linear  combinations  of  u  and  v 


transform  Equations  170  and  171  into 

i(uH-f)v*  -  ((a-b^Ev  ^  +  1(7^^  _  .r8y)j  (172) 

i(u>-f)v-  -  ((a-b^Ey  ^  _  7-^^)  _  1(7.^^  _  7.gy)j  (I73) 


121.  Equations  172  and  173  show  that  the  vertical  structures  of  v*  and  v'  are 
uncoupled  and  that  each  is  analogous  to  the  structure  of  u  in  the  absence  of  the 
Coriolis  force,  except  that  v*  is  forced  by  the  frequency  w  +  /  and  v*  is  forced  by  the 
frequency  w  -  /  Therefore  the  vertical  structures  of  v*  and  v*  will  depend  on  the 
dimensionless  frequencies  fl*  =  0  +  <9^  and  fl*  =  Q  -  respectively,  where 

3^  =  /hVE^.  At  mid-latitudes,  /~  10'<  sec*‘,  giving  the  range  of  10'*  to  10^. 
This  yields  values  for  U*  and  11'  in  the  same  range  as  11;  consequently  the  results 
shown  in  Figures  16-19  are  also  indicative  of  the  performance  of  the  DSS‘  and  the  VS 
when  the  Coriolis  force  is  included  in  the  governing  equations. 

122.  Analytical  solutions  can  be  obtained  for  the  test  problem  for  a  transient 
forcing  by  decomposing  the  forcing  into  its  Fourier  components,  using  the  periodic 
solutions  presented  above  for  each  Fourier  component  and  superimposing  the  resulting 
periodic  solutions.  In  this  section  an  illustrative  set  of  results  for  bottom  stress  are 
presented  for  the  often-used  problem  of  an  instantaneously  imposed  wind  on  an 
initially  quiescent  channel.  Representative  values  of  L  =  100  km,  h  =  50  m, 

Oq  =  0.01,  and  Ezo=  0.5  m*/s  are  used. 

123.  An  instantaneously  imposed  forcing  cannot  be  represented  exactly  by  a 
finite  Fourier  series;  however. 


-liiO-  =  i  +  2  y 

Tg- steady  ^  "  x(  2n-l ) 


(174) 


111 


gives  an  approximation  to  a  square  wave  of  period  T,  as  shown  in  Figure  20.  By 
selecting  T  to  be  larger  than  the  time  required  for  the  basin  to  reach  steady  state 
and  considering  only  the  period  0  <  t/T  <  1,  a  reasonable  representation  of  an 
instantaneously  imposed  wind  can  be  obtained  and  used  to  develop  an  approximate 
analytical  solution.  Sensitivity  analyses  indicated  that  when  50  or  more  terms  were 
used  in  Equation  174,  minimal  change  occurred  in  the  analytical  solution  of  the  basin 
r^ponse  and  any  change  that  did  occur  was  limited  to  times  very  close  to  zero  (i.e., 
on  the  order  of  t/T  <  1  percent).  Seventy-five  terms  (N=74)  were  used  in  Equation 
174  for  the  solution  shown  in  Figure  20  and  the  runs  presented  below. 

124.  The  VS  and  the  DSS'  for  the  transient  test  case  were  obtained  by 
discretizing  Equations  123  and  124  in  time  using  a  Crank-Nicholson  scheme.  As 
discussed  above,  the  analytical  solution  for  U  was  used  to  force  these  equations, 
thereby  eliminating  any  feedback  of  error  from  the  vertical  representation  into  U. 
Figure  21  presents  a  comparison  between  bottom  stresses  obtained  analytically  and 
from  the  VS  and  the  DSS*.  The  DSS*  with  3  LPs  is  quite  close  to  the  analytical 
solution  except  very  near  t  =  0  (due  primarily  to  the  overshoot  in  the  forcing  in 
Figure  20).  Conversely,  15  or  more  LPs  are  required  for  the  VS  to  attain  comparable 
accuracy.  We  note  that  this  test  case  uses  (t©  at  the  upper  limit  of  the  practical 
range  and  therefore  is  the  easiest  case  for  the  VS  to  capture.  For  smaller  values  of 
(To,  the  transient  performance  of  the  VS  becomes  even  poorer  as  suggested  by  the 
steady-state  results  in  Table  4. 

125.  The  results  of  this  section  suggest  that  shear  stress  can  be  a  highly 
efficient  substitute  for  velocity  as  the  dependent  variable  in  the  internal  mode 
equations.  For  this  to  be  accomplished,  it  is  only  necessary  that  the  shear  stress  and 
the  vertical  gradient  of  velocity  be  linked  via  an  eddy  viscosity  relationship. 

Depending  on  the  choice  of  shape  functions  and  the  functional  variation  of  eddy 
viscosity  over  the  depth,  the  velocity  profile  can  be  recovered  from  the  stress  profile 
in  closed  form.  Under  these  conditions  the  difficulties  associated  with  numerically 
integrating  a  near-logarithmic  singularity  are  avoided.  Most  practical  problems  can  be 
solved  subject  to  this  restriction  by  allowing  a  global  or  piecewise  polynomial  variation 
of  Tz  and  a  piecewise  linear  variation  of  Eg. 

126.  One  disadvantage  with  the  DSS'  is  that  it  yields  a  fully  populated  matrix 
on  the  left  side  of  the  discretized  equations  that  must  be  reformed,  decomposed,  and 
solved  at  every  time  step  if  a  time-varying  eddy  viscosity  is  used.  This  requires  ~ 
operations  to  solve  for  stress  and  ~  N*  operations  to  extract  velocity  (using 
Equation  118),  where  N  is  the  number  of  LPs  that  are  used.  Although  often  only  a 
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t  (hrs)  t  (hrs) 

^  nmncrical  and  analytical  bottom  stress  for  the  transient  test  case  for  the 


few  LPs  are  required  for  an  accurate  solution,  as  N  reaches  ~  10,  the  computational 
attractiveness  of  the  DSS‘  rapidly  diminishes  in  comparison  to  a  VS  that  only  requires 
the  solution  of  a  banded  matrix  (e.g.,  Lynch  and  Werner  1991).  Part  of  the  reason 
for  the  fully  populated  matrix  is  due  to  the  spectral  method’s  use  of  globally,  rather 
than  locally,  defined  functions.  If  Equation  124  is  discretized  using  the  finite  element 
method  with  linear  elements,  the  left-side  matrix  is  the  sum  of  a  triangular  plus  a 
tri-diagonal  matrix.  This  requires  ~  operations  to  solve,  where  M  is  the  number 
of  nodes  used  over  the  depth.  It  can  be  shown  that  the  triangular  part  of  the  matrix 
arises  because  of  the  integral  term  in  Equation  124. 

Development  and  Testing  of  DSS  Method  No.  2 


127.  The  solution  of  a  fully  populated  or  near-triangular  matrix  system  can  be 
avoided  by  reformulating  the  DSS  internal  mode  equations  to  eliminate  integral  terms 
£rom  the  left  side.  This  can  be  accomplished  by  generating  internal  mode  equations 
by  taking  the  vertical  derivative  of  the  three-dimensional  momentum  equations  rather 
than  by  subtracting  the  vertically  integrated  equations  from  the  three-dimensional 
equations.  The  use  of  internal  mode  equations  derived  by  taking  the  vertical 
derivative  of  the  three-dimensional  equations  has  been  reported  by  Tee  (1979). 
Although  this  report  focuses  primarily  on  the  simplified  internal  mode  equations  for  a 
constant-density  fluid,  the  derivation  of  the  full  internal  mode  equations  is  presented 
below  for  completeness. 

128.  Differentiating  the  ^-coordinate  horizontal  momentum  equations 
(Equations  19  and  20)  with  respect  to  a,  and  substituting  Equation  21  for  dp/dcr 
gives 


(Note,  this  is  illustrated  for  the  x-momentum  equation  only.  The  y-momentum 
equation  follows  directly.) 

Using  the  eddy  viscosity  relationship  for  r^x  and  rzy  (Equation  34)  the  vertical 
gradient  of  velocity  can  be  expressed  in  terms  of  the  shear  stress  as 

da  _  Etzx 
W  -  EvCa=b)po 
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(176b) 


dv  _  Hrav 
'Ua  iSv(a-b)po 

Substituting  Equation  176  and  the  expansions 

-dBm  -  ^ W 

dfda\  _  d  f„du\  5u  dv 

^  ■  -g?  ^ 

into  Equation  175  gives: 


H'^ZX 
2v(a— b)po 


fHrzv  (a-b)  d^r^x  _  d 

iv(  a-b)po  Hpo  dfT^  m 


ife]  ■ 


vHTzx 
^v(  a— b)po 


_  ^  r  wHtzx  1  .  Htzz  dv 
do j_Ev(  a— b)po J  Ev(  a— b)/Jo 

Using  the  additional  expansions 

^  r  wHTzx  1  _  WH  dfTzx\  _ 

gg[Ev(  a-b)poJ  (a— b)po 

d  r  Htzx  1  _  H  d  /  T zx\  , 

'gt[Ev(a-b)poJ  (a-b)po  ^Ev  '  ^  I 

d  \  uHtzx  1  _  uH  g  ,Tzx\  . 


iL-  ..  _  mi  o,rzx\  ,  Tzx  PUn 

c5c[Ev(  a-bjpoj  (a-b)po  ^E^^  Ev(a-b)po  dx 


Hrzv  gu  gbx  ,  gmx 
[a-b)po  dy  da'  da 


Tzx  \dE  ,  dnU  ,  gvHl 

+  gr~  + 

ZX 

-b)Po  di 


d  vHTzx  1  _  VH  g  ,Tzx\  . 

g5^lE;-(a-b>oJ  “  (a=b7^  o5^Er^  + 


Tzx  dvE 
2v(a-b)po  W~ 


Equation  177  can  be  written  in  final  form  as 


g  Ttzx  1  _  fT'zv  (a-b)^  d^Tzx  _  _  (a — ^b)fgbx  3mxl 

EtJo  ~gg5-  -  cx  -  -  -^J 

where  Cx  represents  the  contribution  of  the  nonlinear  advective  terms 

=  _  ii^FXsilJ  _  \T,t^  1  -  f  '^zx  1  I  Tzx  dv  T  zv  gu 

e;;^  ^  "  e^  ^ 

Applying  the  same  transformation  to  the  y-momentum  equation  gives 
d  T zv  I  f^zx  (a~b)^  d^Tzy  _  (a — b)r^v  gmv 


fa-bir^. 


where 


,,d  \Tzv 


.d  fr 


d  \t 


(177) 


/•  I  ..V  T %y  d  T xy  I  T XV  ^  T xx  dv 

'7  -  -  ® 

Introducing  the  complex  shear  stress  Tz  =  t^x  +  irzy  (where  i  =  ^  -'1  ), 


(178) 


(179) 


(180) 


(181) 
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Equations  178  and  180  can  be  combined  into  a  single  complex  equation 

c  +  b  +  m 


d 


Tz 


(a-b)^  d^Tz 


where 

C  —  Cx  "I"  3iCy 


If  +  i  i"' 


m  = 


ftllx  1  :  ftHv 

3^  +  ‘  3? 


(182) 

{183a) 

(183b) 

(183c) 


Because  both  r®  and  Ey  may  vary  in  time,  the  discretization  of  Equation  182  in  time 
may  be  facilitated  by  expanding  the  leading  term  as: 


d  Tz  _  1  d  fTz\  _  Tz  1  ^Ey 

^[EyPoJ  ~  11^  '^po'  Po'E^  M 


(184) 


Substituting  Equation  184  into  Equation  182  and  multiplying  both  sides  by  Ey  gives 


5Ey 

'diSpo'  Po^y  ^ 


+  if^  - 

Po  B  ^Pq 


-  Ev[c  +  6  +  m] 


(185) 


129.  For  the  3DL  model,  the  baroclinic,  advective,  and  horizontal  turbulent 
momentum  terms  are  assumed  to  be  equal  to  zero.  This  leaves 


d 

'di 


4.  if 


(a-b)2  d^Tz 


0 


(186a) 


(186b) 


as  simplified  DSS^  internal  mode  equations.  (The  superscript  2  is  used  to  identify 
DSS  method  No.  2.)  We  note  that  for  an  eddy  viscosity  that  is  constant  in  time, 
Equations  186a  and  186b  have  the  form  of  complex  diffusion  equations  for  stress. 

This  provides  a  physical  interpretation  for  the  internal  mode  equation;  i.e,  it  describes 
the  turbulent  diffusion  of  stress  through  the  water  column. 

130.  Because  of  the  second  derivative  term  in  stress  in  Equations  186a  and 
186b,  two  boundary  conditions  are  needed  to  solve  either  equation  over  the  vertical. 
The  free  surface  boundary  condition  is 


Tz/po  =  Ts/po  at  <r  =  a 


(187) 


where  Tg  is  the  specified  surface  stress.  A  second  boundary  condition  can  be 
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generated  by  requiring  the  depth  average  of  the  internal  mode  velocity  to  match  the 
external  mode  depth-averaged  velocity.  From  Equation  118  this  becomes 

+  (5^.  ITSe;  drf.  =  u  +  iv  (188) 

b  b 


131.  To  avoid  the  fully  populated  matrices  generated  by  the  Galerkin-spectral 
method,  the  DSS®  uses  the  Galerkin-finite  element  method  to  discretize  the  internal 
mode  equation  over  the  vertical,  r^/po  is  expanded  over  M-1  depth  intervals  using 
depth-dependent,  locally  defined  basis  functions  Fr(<r)  and  complex  coefficients  7r(x,y,t) 

u 

Po  “ 


The  Galerkin-finite  element  forms  of  Equations  186a  and  186b  are  obtained  by 
substituting  Equation  189  for  Tz/po,  multiplying  each  equation  by  Fr(<T)  and 
integrating  with  respect  to  a  over  the  interval  from  a  to  b: 

^  f  a  a  *  .a 


F  dn 


m*l  L 

and 

u 


=  0 


n  =  1,  ...M 


I  I?  +  -  7i»  I  ^  da  -  7„ 

■  =  iL  i  i  '' 


F  F  dn 


n  =  1,  ...M 


(190a) 

=  0 
(190b) 


132.  Linear  chapeau  functions  will  be  used  for  Fr(a).  The  tendency  observed 
in  the  DSS*  results  for  stress  to  become  linear  over  the  depth  for  fi  <  1  suggests  that 
these  functions  should  give  a  good  representation  of  stress  if  the  element  size  is 
selected  so  that  n®  <v  1.  (Qe  is  identical  to  Q  except  it  is  defined  using  the  element 
size  rather  than  the  total  depth.)  However,  Equations  190a  and  190b  require  a 
interpolating  basis.  To  lower  this  requirement  to  C°,  we  integrate  by  parts: 


I 

I 


Fn 


da  = 


i 

-I 


dFn  dFn 

W  W 


da 


(191a) 


118 


(191b) 


j  E.fn  <i»  =  Ev(a)f„(a)5^  -  E.(b)F„(b)^^  -  |  |5(EvF.) 


for  Equations  190a  and  190b,  respectively. 

133.  Using  linear  basis  functions,  when  n  =  1  and  n  =  M,  the  first  two  terms 
in  Equations  191a  and  191b  exactly  cancel  the  integral  terms  in  these  equations 
making  the  total  diffusion  terms  equal  to  zero.  However,  when  2  <  n  <  M-1,  the 
first  two  terms  in  Equations  191a  and  191b  are  identically  zero.  Therefore,  for 
2  <  n  <  M-1,  Equations  191a  and  191b  can  be  substituted  into  Equations  190a  and 
190b  to  give  physically  meaningful  equations: 

V  \d\^  rVaFn  .  -  (a-b)2r^dFn 

Z  ■^t  [^“  J  'ETT  ^  J  J  ^  \  'ST '5a~  ^ 

^  b  b  b  '' 

n  =  2,  ...M-1  (192a) 


"  [■  *  *  2  a 

I  [^t^  ^  -  7ai^|  ^oi-Pnln(Ev)  da  +  7B^^ga^  da  =  0 


n  =  2,  ...M-1  (192b) 

134.  The  boundary  conditions  are  used  to  supplement  Equations  192a  and  192b 
when  n  =  1  or  n  =  M.  Equation  187  is  used  in  place  of  the  n  =  M  equation: 


Re{7y}  =  ’■sx/Po  and  lB{7y}  =  rgy/po 
In  place  of  the  n  =  1  equation,  Equation  188  gives 

^f2iE|Ibl  +  I  da  dal  =  U  +  iV 


(193) 


(194) 


135.  Velocity  is  recovered  from  stress  by  solving  the  discretized  version  of 
Equation  118 


»  +  iv  =  I  2!!E|M  +  I  I® 


(195) 


136.  Equations  192a  or  192b  and  193  form  a  tri-diagonal  system;  Equation  194 
adds  a  fully  populated  bottom  row  to  this  system.  However,  only  a  few  extra 


computations  are  required  to  reduce  the  system  to  tri-diagonal.  Therefore,  the 
number  of  operations  required  to  obtain  a  solution  for  stress  at  each  time  step  scales 
with  M.  Since  T^fpo  is  piece  wise  linear  over  the  depth,  the  integrals  in  Equation 
195  can  be  evaluated  analytically  for  many  functional  forms  of  Ey.  For  most  practical 
model  applications,  it  can  be  assumed  that  Ey  has  a  piece  wise  linear  variation  with 
depth  (Fumes  1983;  Chu,  Liou,  and  Flenniken  1989;  Jenter  and  Madsen  1989)  This 
is  physically  correct  near  boundaries  and  makes  the  analytical  solution  of  the  stress 
integrals  particularly  simple.  Using  this  functional  form  for  Ey,  the  number  of 
operations  required  to  analytically  extract  velocity  from  stress  also  scales  with  M. 

137.  An  initial  evaluation  of  the  DSS^  has  been  made  using  the  same  test 
problems  solved  for  the  DSSi  and  the  VS.  For  these  tests  f  =  0,  Ty=  0,  and  Ey  is 
constant  in  time.  Therefore,  Equations  192  -  195  are  simplified  to: 


1  ^ 


=  0  n  =  2,  ...M-1 


r  I?  I  j  d(7 


=  0  n  =  2,  ...M-1 


7y  —  TsxIPo 


n  =  M 


I  7m 


Fjb)  .  h  f  f 
?ok  (a-b)2j  J  E7 

b  b 


da  d(T 


=  U 


n  =  1 


u  =  Jtb 


n-l 


Pak  ^  (a=BJ  J  e:: 


da 


(196a) 

(196b) 

(197) 

(198) 

(199) 


138.  To  distinguish  between  the  two  internal  mode  equations,  results  are 
designated  as  DSS|  or  DSS|  depending  on  whether  they  are  based  on  Equation  196a 
or  Equation  196b,  respectively.  In  all  of  the  results,  a  specified  number  of  equal-sized 
elements  was  used  over  the  vertical.  It  may  be  possible  to  improve  the  efficiency  of 
the  DSS2  further  using  elements  that  are  not  equally  sized.  However,  this  option  has 
not  yet  been  investigated  completely. 

139.  In  the  steady-state  test  case,  the  stress  distribution  is  linear  over  the 
depth  (Equation  143);  therefore,  both  the  DSSi  and  the  DSS§  give  the  exact  solution 
using  one  element  over  the  vertical.  The  number  of  degrees  of  freedom  (NDF)  in  the 
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finite  element  solution,  (i.e.,  the  number  of  simultaneous  equations  that  must  be 
solved)  is  equal  to  the  number  of  nodes  used  in  the  discretization  (number  of  nodes  = 
1  +  niimber  of  elements).  The  NDF  in  the  spectral  solution  is  equal  to  the  number 
of  LPs  used  in  the  discretization.  In  both  the  DSS^  method  and  the  DSS^  method, 
the  exact  steady-€tate  solution  is  obtained  using  two  degrees  of  freedom. 

140.  Results  from  the  periodic  test  case  are  shown  in  Figures  22  -  26.  When 
ft  <  1,  the  DSS*  is  nearly  exact  using  one  finite  element  (two  degrees  of  freedom) 
over  the  depth,  Figures  22a  and  22b,  23a  and  23b,  24a  and  24b,  and  25a  and  25b. 

For  ft  >  1,  more  than  one  finite  element  is  required  over  the  vertical  for  either  DSS^ 
to  converge  to  the  analytical  solution  (Figures  22c  and  22d,  23c  and  23d,  24c  and 
24d,  and  25c  and  25d).  Comparing  these  results  to  the  DSS<  results  indicates  that 
both  DSS^  methods  require  more  degrees  of  freedom  than  the  DSS»  method  to  reach 
the  same  level  of  convergence.  The  bottom  stress  plots  presented  in  Figure  26 
demonstrate  the  properties  of  the  DSS*  method  further.  In  particular,  they  indicate 
that  the  DSS^  is  quite  effective  in  the  range  ft  ~  10  or  less.  It  may  be  possible  to 
extend  this  range  to  higher  values  of  ft  if  an  unequally  spaced  finite  element  grid  is 
used  over  the  depth. 

141.  A  time  history  of  bottom  stress  for  the  transient  test  case  is  shown  in 
Figure  27.  Comparing  this  to  Figure  21a  indicates  that  both  DSS*  methods  require 
four  degrees  of  freedom  to  give  a  solution  that  is  approximately  equivalent  to  the 
DSS‘  using  three  degrees  of  freedom. 

142.  In  conclusion,  new  internal  mode  equations  have  been  developed  that 
allow  shear  stress  to  be  used  as  the  dependent  variable  in  the  internal  mode  solution 
and  that  yield  a  nearly  tri-diagonal  matrix  system.  While  both  DSS*  require  more 
degrees  of  freedom  than  the  DSS^  method  to  obtain  comparable  results  for  ft  >  1, 

(due  to  the  use  of  linear  finite  elements  in  the  DSS*  versus  spectral  functions  in  the 
DSSi),  the  matrix  structure  of  the  DSS*  matrices  makes  this  method  much  more 
efficient  than  the  DSS‘. 

Implementation  of  Wave-Current  Interaction  in  a  DSS  Model 

143.  It  is  often  observed  in  lakes,  coastal  waters,  and  shelf  waters  that  near 
the  bottom  the  orbital  velocities  associated  with  surface  waves  are  as  large  as  or 
larger  than  the  mean  current  velocity.  In  such  cases  the  surface  waves  have  a 
significant  effect  on  the  bottom  stress  and  the  current  profile.  Several  investigators 
have  developed  theoretical  models  to  account  for  this  wave-current  interaction.  To 
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Figure  24.  Vertical  profiles  of  horizontal  velocity  magnitude  for  the  periodic  test  case  using  the  DSSi 


Vertical  profiles  of  horizontal  velocity  phase  for  the  periodic  test  case  using  the  DSS^ 
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thoroughly  assess  the  usefulness  of  the  DSS  approach,  the  effort  required  to  implement 
the  Grant  and  Madsen  (1979)  model  (GM  model)  with  a  DSS  of  the  internal  mode 
equations  has  been  considered.  The  GM  model  assumes  that  the  mean  current 
velocity  can  be  determined  as  follows: 

a.  inside  the  wave  boundary  layer,  z  <  6m,, 

Ey  =  K  jU^icwl  2  (200) 

ju+cwl  =  ^  (201) 

a  no-slip  boundary  condition  is  applied  at  z  =  Zo,  where  Zq  is  the 
physical  bottom  roughness 

b.  outside  the  wave  boundary  layer,  z  >  S,, 

Ey  =  «  |U*c|  z  (202) 

iu.c|=ijrcr=  (203) 

a  no-slip  boundary  condition  is  applied  at  z  =  zoa,  where  Zoa  is  an 
apparent  bottom  roughness  experienced  by  the  current  due  to  the 
wave-current  interaction. 

In  these  relations,  ^  =  0.4  is  the  Von  Kirman  constant,  Tc  is  the  bottom  stress  due 
to  the  current  alone,  Tw  is  the  maximum  wave-4nduced  bottom  stress  during  a  wave 
cycle,  and  is  the  thickness  of  the  wave  boundary  layer. 

144.  The  GM  model  can  be  included  in  a  DSS  of  the  internal  mode  equations 
as  follows. 

g.  Estimate  Zoa  and  jU^d  based  on  values  at  the  previous  time  step. 

b.  Calculate  Ey  and  use  the  DSS  model  to  predict  Tc- 

c.  Solve  Equation  201  for  lU*cw!  using  Tq  from  the  previous  step  and 

Tm,  from  Equation  53  in  Grant  and  Madsen  (1979).  Since  Tw  is  a 
function  of  Equation  201  must  be  solved  iteratively. 

d.  Determine  Zoa  using  Equations  46  and  49  in  Grant  and  Madsen 
(1979). 

g.  Recalculate  Ey  using  the  new  Tq.  Use  tliis  and  the  new  value  of  Zoa 
in  the  DSS  model  to  predict  Tc-  Go  to  step  c.  and  iterate  until  Tq 
(X)nverges. 

145.  Because  two  levels  of  iterations  are  required  to  implement  the  wave- 
current  interaction,  it  may  be  computationally  infeasible  to  use  this  scheme  in 
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practical  model  applications.  It  may  be  possible  to  simplify  this  procedure  in  two 
ways  in  the  proposed  model.  First,  rather  than  iterate  as  described  in  step  g.,  Zoa 
can  be  calculated  explicitly  in  time  based  only  on  results  from  the  previous  time  step. 
This  should  introduce  little  error  into  the  solution  if  the  time  step  is  small  enough 
that  changes  in  Zoa  and  Tq  are  relatively  small.  Second,  following  the  suggestion  of 
Spaulding  and  Isaji  (1987),  Tw  can  be  determined  by  neglecting  the  effect  of  the 
current  on  the  wave  within  the  wave  boundary  layer.  In  this  case 

Tw/p  =  0.5  fw  lUbI  Ub  (204) 

where  fw  is  the  wave  friction  factor  (Jonsson  and  Carlsen  1976)  and  Ub  is  the 
maximum  bottom  wave  orbital  velocity.  For  fully  rough,  turbulent  flow,  fw  can  be 
determined  from 

— 1=  +  l08.»  =  ‘°6io  ^  -  0.12  (205) 

4  ^  fw  4  ^  fw 

where  Ab  is  the  bottom  excursion  amplitude  of  the  wave  and  kg  is  the  Nikuradse 
equivalent  sand  roughness  of  the  bottom  (typically  Zq  =  ks/30). 

146.  The  brief  outline  presented  above  suggests  that  the  GM  wave-current 
interaction  can  easily  be  included  in  the  DSS  model.  In  fact,  if  the  implementation 
procedure  outlined  above  for  the  DSS  is  compared  with  that  described  in  Grant  and 
Madsen  (1979)  for  a  standard  VS,  it  is  evident  that  the  DSS  simplifies  the  use  of  the 
GM  model  by  eliminating  the  complications  introduced  by  a  quadratic  slip  bottom 
boundary  condition. 


129 


PART  V:  SUMMARY  AND  CONCLUSIONS 


147.  This  report  documents  the  theory  and  methodology  behind  the  ADCIRC 
(Advanced  Circulation)  model’s  2DDI  (2-dimensional,  depth-integrated)  option  and  the 
3DL  (3-dimensional,  local  internal  mode  equation)  option.  ADCIRC  is  based  on  the 
three-dimensional  Reynold’s  equations  simplified  using  the  hydrostatic  pressure  and  the 
Boussinesq  approximations.  Prior  to  their  solution,  the  three-dimensional  equations 
are  separated  into  a  set  of  external  mode  equations  (the  two-dimensional,  vertically 
integrated  equations)  and  a  set  of  internal  mode  equations. 

148.  The  external  mode  equations  can  be  solved  by  themselves  (the  2DDI 
option)  for  depth-averaged  velocity  and  free-surface  elevation  by  parameterizing 
bottom  stress  and  momentum  dispersion  in  terms  of  the  depth-averaged  velocity.  Key 
features  of  the  external  mode  solution  are  the  use  of  a  generalized  wave-continuity 
equation  (GWCE)  formulation  and  the  Galerkin-finite  element  (FE)  method  in  space 
using  triangular  or  quadrilateral  elements.  The  FE  method  provides  maximum  grid 
flexibility  and  allows  highly  efficient  numerical  solutions  to  be  obtained  using  model 
domains  that  include  complicated  bathymetries  and  shoreline  geometries  that  also 
stretch  considerable  distances  offshore  to  implement  open-water  boundary  conditions. 
Detailed  analyses  and  testing  of  ADCIRC-2DDI  have  shown  that  it  has  good  stability 
characteristics,  generates  no  spurious  artificial  modes,  has  minimal  inherent  numerical 
damping,  and  efficiently  separates  the  external  mode  equations  into  small  systems  of 
algebraic  equations  with  time-4ndependent  matrices.  Applications  of  the 
ADCIRC-2DDI  model  to  the  English  Channel  and  southern  North  Sea,  the  Gulf  of 
Mexico,  Masonboro  Inlet,  and  the  New  York  Bight  have  shown  that  it  is  capable  of 
running  month  to  year-long  simulations  while  providing  detailed  intra-tidal 
computations. 

149.  In  stratified  flows,  Ekman  layers,  wind-driven  flows  in  enclosed  or  semi- 
enclosed  basins,  or  flows  affected  by  wave-current  interaction  in  the  boundary  layer,  it 
is  generally  impossible  to  parameterize  bottom  stress  and  momentum  dispersion  in 
terms  of  depth-averaged  velocity.  In  such  cases,  it  is  necessary  to  solve  the  internal 
mode  equations  for  the  vertical  variation  of  horizontal  velocity  and  use  this  to 
evaluate  the  bottom  stress  and  momentum  dispersion  terms  in  the  external  mode 
equation.  Due  to  the  shallow  water  depths  that  characterize  coastal  and  shelf 
settings,  the  internal  mode  equations  can  often  be  simplified  by  dropping  the 
horizontal  gradient  terms.  This  gives  internal  mode  equations  that  express  the  vertical 
distribution  of  momentum  at  any  horizontal  position  as  a  local  balance  between  the 
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surface  and  bottom  stresses,  vertical  momentum  diffusion,  the  Coriolis  ^orce,  and  local 
inertia.  The  SDL  model  option  is  formulated  using  the  simplified,  local  internal  mode 
equations.  Existing  numerical  solutions  of  full  or  simplified  internal  mode  equations 
use  velocity  as  the  dependent  variable.  Therefore,  it  is  necessary  to  use  a  fine 
numerical  discretization  to  resolve  the  sharp  vertical  gradients  of  velocity  that  occur 
near  the  bottom  boimdary  and  in  wind-driven  flows  near  the  surface  boundary. 

During  the  course  of  the  ADCIRC-3DL  model  development,  a  novel  technique  was 
discovered  that  replaces  velocity  with  shear  stress  as  the  dependent  variable  in  the 
internal  mode  equations.  The  resulting  direct  stress  solution  (DSS)  allows  physically 
realistic  boundary  layers  to  be  explicitly  included  in  a  three-dimensional  model. 
Detailed  testing  of  the  DSS  method  has  demonstrated  its  considerable  advantage  over 
standard  velocity  solutions  and  has  led  to  an  optimized  DSS  formulation.  This 
treatment  of  the  internal  mode  equations  should  be  invaluable  for  modeling  coastal 
and  shelf  circulation  in  which  the  bottom  and  surface  boundary  layers  comprise  a 
significant  portion  of  the  water  column  and  for  modeling  processes  that  are  critically 
dependent  on  boundary  layer  physics  such  as  wave-current  interaction,  sediment 
transport,  oil  spill  movement,  ice  floe  movement,  energy  dissipation,  physical-biological 
couplings,  etc. 

150.  Considerable  effort  has  gone  into  the  development  of  ADCIRC  to  produce 
a  model  that  has  simultaneous  regional/local  capabilities,  as  well  as  very  high  levels  of 
accuracy  and  efficiency.  This  has  been  achieved  by  combining  extreme  grid  flexibility 
with  optimized  formulations  of  the  governing  equations  and  numerical  algorithms. 
Together,  these  allow  ADCIRC  to  run  with  improved  physical  realism  and  a 
significant  reduction  in  the  computational  cost  of  most  presently  existing  circulation 
models. 
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